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PREFACE 


This volume of the Annual Status Report describes the results of work 
performed on Task IVB (Special Finite Element Models) during the third year of 
the NASA Hot Section Technology Program, "3-D Inelastic Analysis Methods for 
Hot Section Components" (Contract NAS3-23697). The goal of the program is to 
develop computer codes which permit more accurate and efficient structural 
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PROLOGUE 


3-D INELASTIC ANALYSIS METHODS FOR HOT SECTION COMPONENTS 
VOLUME I - SPECIAL FINITE ELEMENT METHODS 


by 

S. Nakazawa 

MARC Analysis Research Corporation 
Suite 314 Court House Plaza, 260 Sheridan Avenue, 
Palo Alto, California 94306 


The Special Finite Element Models portion of the 3-D Inelastic Analysis 
Methods Program is divided into two 24-month segments: a base program, and an 
option program exercised at the discretion of the Government. Versions 1 and 2 
of the MHOST (MARC - HOST ) code were developed during the base program (Tasks IB 
and I IB ) . The MHOST code employs both shell and solid (brick) elements in an 

iterative mixed method framework to provide comprehensive capabilities for 

investigating local (stress/strain) and global (vibration/buckling) behavior 
of hot section components. In the development of the code, full advantage has 
been taken of the wealth of technical expertise accumulated at the MARC 
Corporation over the last decade in support of their own commercially 

available software package. This has led to the creation of new/improved 

algorithms that promise to significantly reduce the central processing unit 
(CPU) time requirements for three-dimensional analyses. 

Version 3 of the MHOST code concentrated on algorithm refinement to boost the 
performance of the basic mixed iterative solution methods; work has also 
continued on a subelement solution strategy designed to capture local behavior 
related to embedded discontinuities. Considerable efforts have been devoted to 
maintaining the computer program and improving the quality of the code, 
including a substantial number of comment statements in the code. An 
out-of-core frontal solution subsystem has been made available, in addition to 
the original band matrix solution, and the MHOST user interface has been 
modified to make it more user-friendly. An anisotropic material data reader 
has been added to MHOST as a replacement for the user subroutine option. The 
shell element definition has also been modified to allow the modeling of 
composite laminate materials. 
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SECTION 1.0 
INTRODUCTION 


Aircraft powerplant fuel consumption and expenditures for repair/replacement 
of worn or damaged parts make up a significant portion of commercial 
aviation's direct operating costs. For modern gas turbines, both factors 
depend heavily on the degree to which elevated flowpath temperatures are 
sustained in the hot section modules of the engine. Higher temperatures reduce 
fuel consumption by raising the basic efficiency of the gas generator 
thermodynamic cycle. At the same time, these elevated temperatures work to 
degrade the durability of structural components (combustor liners, turbine 
blades and vanes, airseals, etc.) that must function adjacent to or within the 
hot gaspath itself, leading in turn to larger maintenance/material costs. 
Pursuit of the best compromise between performance and durability presents a 
challenge that will continue to tax the ingenuity of advanced gas turbine 
design analysts for years to come. 

Hot section durability problems appear in a variety of forms, ranging from 
oxidation/corrosion, erosion, and distortion (creep deformations) to 
occurrence of fatigue cracking. Even modest changes in shape, from erosion or 
distortion of airfoils for example, can lead to measurable performance 
deterioration that must be accurately predicted during propulsion system 
design to ensure that long-term efficiency guarantees can be met. Larger 
distortions introduce serious problems such as hot spots and profile shifts 
resulting from diversion of cooling air, high vibratory stresses associated 
with loose turbine blade shrouds, difficult disassembly/reassembly of mating 
parts at overhaul, etc. These problems must be considered and efforts made to 
eliminate their effects during the engine design/development process. 
Initiation and propagation of fatigue cracks represents a direct threat to 
component structural integrity and must be thoroughly understood and 
accurately predicted to ensure continued safe and efficient engine operation. 

Accurate prediction of component fatigue lives is strongly dependent on the 
success with which inelastic stress/strain states in the vicinity of holes, 
fillets, welds, and other discontinuities can be calculated. Stress/strain 
computations for hot section components are made particularly difficult by two 
factors: the high degree of geometrical irregularity which accompanies 
sophisticated cooling schemes, and complex nonlinear material behavior 
associated with high temperature creep/plasticity effects. Since cooling air 
extraction reduces engine cycle efficiency, concerted efforts are made to 
minimize its use with the result that elaborate internal passages and surface 
ports are employed to selectively bathe local regions (airfoil leading edges, 
louver liner lips, etc.) for which the high temperature environment is most 
severe. These cooling features frequently interrupt load paths and introduce 
complex temperature gradients to the extent that the basic assumptions of one- 
and two-dimensional stress analysis procedures are seriously compromised and 
the use of three-dimensional techniques becomes mandatory. Even in the 
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presence of cooling, component temperature and stress levels remain high 
relative to the material's melting point and yield strength values. The 
combinations of centrifugal, aerodynamic, thermal, and other mechanical 
loadings that typically occur in flight operation then serve to drive the 
underlying material response beyond accepted limits for linear elastic 
behavior and into the regime characterized by inelastic, time-dependent 
structural deformations. Thus, an ability to account for both complexities, 
three-dimensional and inelastic effects, becomes essential to the desiqn of 
durable hot section components. 

General purpose finite element computer codes containing a variety of 
three-dimensional (brick) elements and inelastic material models have been 
available for more than a decade. Incorporation of such codes into the hot 
section design process has been severely limited by high costs associated with 
the extensive labor/computer/time resources required to obtain reasonably 
detailed results. Geometric modeling systems and automated input/output data 
processing packages have received first attention from software developers in 
recent years and will soon mature to the point that previous over-riding 
manpower concerns will be alleviated. Prohibitive amounts of central 
processing unit (CPU) time are still required for execution of even 
modest-size three-dimension inelastic stress analyses, however, and is chief 
among the obstacles remaining to be remedied. With today's computers and 
solution algorithms, models described by a few hundred displacement 
degrees-of-freedom commonly consume one to three hours of mainframe CPU time 
during simulation of a single thermomechanical loading cycle. A sequence of 
many such cycles may, of course, be needed to reach the stabilized conditions 
of interest. Since accurate idealizations of components with only a few 
geometrical discontinuities can easily contain several thousand 
degrees-of-freedom, inelastic analysis of hot section hardware with existinq 
codes falls outside the realm of practicality. 

The Inelastic Methods Program addresses the need to develop more efficient and 
accurate three-dimensional inelastic structural analysis procedures for gas 
turbine hot section components. A series of new, increasingly rigorous, 
stand-alone computer codes is being created for the comprehensive numerical 
analysis of combustor liners, turbine blades and vanes. One of these 
stand-alone computer programs is based on special finite element models. Heavy 
attention is being given to the evolutionary development of novel finite 
element modeling methods that permit non-burdensome yet accurate 
representations of geometrical discontinuities such as cooling holes and 
coating cracks. A selection of constitutive relations is being provided for 
either economical or sophisticated descriptions of inelastic material behavior 
as desired. The advantages which accrue from application of the improved 
finite element code being developed under this contract are demonstrated by 
execution of benchmark analyses for which experimental data exist as well as 
by efficiency comparisons with other commonly available finite element 
programs. 
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SECTION 2.0 


SUMMARY 


The Special Finite Element Models portion of the 3-D Inelastic Analysis 
Methods program is divided into two 24-roonth segments: a base program, and an 
option program to be exercised at the discretion of the Government. Versions 1 
and 2 of the MHOST (MARC - HOST ) code were developed during the base program 
(Tasks IB and I IB) . The MHOsT code employs both shell and solid (brick) 
elements in an iterative mixed method framework to provide comprehensive 
capabilities for investigating local (stress/strain) and global (vibration, 
buckling) behavior of hot section components. In the development of the code, 
full advantage has been taken of the wealth of technical expertise accumulated 
at the MARC Corporation over the last decade in support of their own 
commercially available software package. This has led to the creation of 
new/improved algorithms that promise to significantly reduce the central 
processing unit (CPU) time requirements for three-dimensional analyses. 

Version 3 of the MHOST code has been developed under Task IVB of the 3-D 
Inelastic Analysis Methods option program. The development efforts have 
concentrated on algorithm refinement to boost the performance of the basic 
mixed iterative solution method; work has also continued on a subelement 
solution strategy designed to capture local behavior related to embedded 
discontinuities. 

Considerable efforts have been devoted to maintaining the computer program and 
improving the quality of the code, including a substantial number of comment 
statements in the code. An out-of-core frontal solution subsystem has been 
made available, in addition to the original band matrix solution, and the 
MHOST user interface has been modified to make it more user-friendly. 

At the request of the NASA Program Manager, an anisotropic material data 
reader has been added to MHOST as a replacement for the user subroutine 
option. The shell element definition has also been modified to allow the 
modeling of composite laminate materials. 

Algorithmic development for the mixed global solution strategy and local 
solution techniques as well as computer program development are discussed in 
detail in Section 3 of this report. 
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SECTION 3.0 


TECHNICAL PROGRESS 


3.1 Literature Survey 

This section reviews currently available literature on numerical technology 
with special emphasis on the application of finite elements to nonlinear 
problems. The main objective is to survey the present state of the art. 

Papers, reports and books published after 1980 are included in this survey in 
conjunction with a few historically significant references published earlier. 

3.1.1 Introduction 

The finite element method has become a standard tool for numerically 
simulating a wide range of engineering problems using well-established 
technology and commercially available software. The numerical methods have 
been developed continuously in parallel with progress in mathematical theory 
and hardware technology. The solution of three-dimensional nonlinear problems 
still remains costly. Analyses with sufficient resolution could be 
prohibitively expensive without major improvements in numerical and 
computational methodology. 

A major thrust in research and development in the present decade is to refine 
the algorithms and approximations in order to improve the computational 
performance of the finite element method. A significant development in the 
last decade is the penetration of this numerical technology into new 
engineering fields of applications as summarized in the treatise by 
Zienkiewicz (1977). This book is currently being revised to include some of 
the major results of the present decade. 

Problems tackled by the finite element method which require new developments 
are generally nonlinear, multi-dimensional and often subjected to complicated 
constraint equations and inequalities. Standard technology is often 
inappropriate or impractical due to the limitations of hardware performance: 
speed and storage capacity. Research efforts concentrated on the following 
areas to fully exploit the computational resources: 

1. Variational formulations and discretization processes 

Organized "variational crimes" were committed in the past just to 
make the displacement method operational. These crimes are now 
legalized by virtue of various mixed and hybrid formulations. The use 
of the Hu-Washizu principle now provides a firm variational 
foundation for modern finite element technology. The patch test has 
also been investigated in the mixed and hybrid methods contexts in 
recent literature. Augmented forms of the variational statement, 
which are used to construct numerical solution processes for 
constrained problems, appear in modern finite element literature. 
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2. Global solution procedures for linear and nonlinear finite element 
equations 

As an alternative to factorization of the global stiffness equations, 
the utilization of iterative solution algorithms has gained 
popularity. The research effort is concentrated on constructing 
inexpensive, yet accurate preconditioners in such processes. The 
iterative solution methods developed for linear and nonlinear 
programming are used to solve nonlinear finite element equations, as 
found in papers and reports published recently. 

3. Postprocessing techniques including stress smoothing, a posteriori 
error estimates and adaptive mesh refinement 

Technology to extract extra information from a finite element 
displacement solution has been developed with a focus on algorithms 
designed to recover the stress accurately. In addition, research in 
the underlying variational theory has made notable progress. Error 
estimates and indicators based on the solution are now available with 
rigorous mathematical theory. These postprocessing techniques lead to 
the development of mechanisms to adaptively refine finite element 
meshes. 

Branches of engineering in which the finite element method is applied have 
diversified to include fluid mechanics and transport phenomena of heat and 
mass, among many others. The emphasis in this report is on solid and 
structural mechanics including basic mathematical and numerical ideas and 
technology specific to this field. Application to three-dimensional inelastic 
solutions for turbine engine hot section components subjected to realistic 
thermal and mechanical loading cycles is of special interest. The main cause 
of difficulties here involves the complicated material response under high 
temperature and pressure. As summarized by Thompson (1982), the constitutive 
laws built for monolithic alloys are not simple under practical operating 
conditions. Introduction of composites such as fiber reinforced superalloys 
adds further complications to the numerical modeling procedure [Hopkins and 
Chamis (1984)] and requires an enhancement of present methods to handle such 
material models. The physical background of the problem is summarized in 
Hartman (1984), and Cook, Laflen (1984). 

A literature survey by Fyhre and Hughes (1983) covers the papers and reports 
on finite elements potentially useful to this class of problems. The report 
presents a broad overview of the field rather than indicating a specific 
direction. The objective of this survey is more specific. After the 
development of new technology reported in Todd, Cassenti, Nakazawa and 
Banerjee (1984), efforts are now concentrated on the numerical techniques 
directly applicable to further enhancement of the methodology generated in 
this project. 
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3.1.2 Finite Element Approximations 

This section surveys the development of spatial discretization nethodoloq.y in 
finite elements. Historically, numerical tricks and local approximations have 
been introduced in the displacement finite element method. These improve the 
numerical performance, but do not always satisfy the variational principle. 
Thus, "variational crimes" originated. Then, efforts to mathematically 
legitimize such ad hoc schemes followed. 

The variational formulation and the element technology are not clearly 
separable. However, to make the underlying mathematical structure visible, the 
following sections attempt to separately review: (1) the global variational 
forms, and (2) the construction of elements. 

3. 1.2.1 Variational Formulations 

The most general form of the variational functional was constructed oriqinally 
by Hu (1955) and Washizu (1955) for a class of problems in linear elasticity. 
The development of this principle for a range of problems in mechanics is 
fully documented by Washizu (1974) and Oden, Reddy (1976), and for the 
derivation for linear elasticity by Oden (1979). The practical importance of 
this form was not recognized until recently [Zienkiewicz (1977), ch. 12, page 
310]. 

The historical development of variational calculus and the important 
fundamental theorems are concisely presented in Goldstein (1980). For 
application to the finite element method, the conceptual background can be 
traced back to the paper presented by Courant (1943). 

Recent publications address the importance of the Hu-Washizu form and its 
applications. Zienkiewicz and Nakazawa (1984) discuss the equivalence of the 
Hu-Washizu principle and numerically integrated displacement finite element 
formulations. This is an important generalization of the equivalence theorem 
by Malkus and Hughes (1978), originally stated for a simpler class of 
constrained problems. Hughes and Malkus (1983) present another extension of 
the theorem, still in the framework of the Hel linger-Reissner principle, to a 
class of anisotropic and nonlinear problems. 

Loubignac, Cantin, Touzot (1977) and Cantin, Loubignac, Touzot (1978) 
introduce a method to solve elastic problems iteratively by using continuous 
stress approximations. The procedure is a different manifestation of the mixed 
method derived from the Hel linger-Reissner principle [Zienkiewicz, Li and 
Nakazawa (1984)]. This procedure indeed improves the quality of the the finite 
element solution, in particular the stress field. A more accurate version is 
reported by Zienkiewicz, Vilotte, Toyoshima and Nakazawa (1984). It compares 
various integration schemes for the constraint equation and finds that some of 
them produce simultaneously super-convergent displacement and stress at nodes. 
The paper also points out the similarity of this procedure and the augmented 
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Lagrangian iteration. An extension of this procedure for a general class of 
inelastic problems requires the full Hu-Washizu principle. A first draft of 
this discussion has appeared in Todd, Cassenti, Nakazawa and Banerjee (1984) 
as a contract report. Construction and utilization of the augmented Hu-Washizu 
form and iterative solution algorithms are discussed further in Nakazawa, 
Nagtegaal and Zienkiewicz (1985). 

An iterative cross section method proposed by Crisfield (1975) looks familiar 
to the algorithm discussed above. The motivation was to improve in an 
iterative manner the transverse shear constraint for the incompatible plate 
element iteratively. As later criticized by Mang and Gallagher (1978), the 
variational structure of this scheme is not consistent with the theory and the 
result is dubious. The construction of iterative solution algorithms for 
incompatible plates with a firm variational foundation first appeared in Todd, 
Cassenti, Nakazawa and Banerjee (1984). In addition. Miller and Hughes (1985) 
discuss a version of a consistent algorithm for structures. 

Development of a generalized nonlinear version of the Hu-Washizu principle for 
finite strain, incompressible plasticity is reported by Simo, Taylor and 
Pister (1984). The principle is also extended to the formulation of finite 
strain viscoelasticity by Simo (1985). 

The Hu-Washizu principle is used extensively as a vehicle to derive the hybrid 
variational principle and the hybrid finite elements. An outline of the 
derivation is presented in Oden and Reddy (1976). An example of the hybrid 
plate stiffness is given in Spilker and Belytschko (1983). The finite element 
formulation is similar (often identical) to the discretized plate and shell 
model originally derived by Reissner (1948) and Mind 1 i n (1951). These forms 
relax the C^ -continuity requirement of the admissible variations and make 
the construction of finite elements simpler. The numerical methods derived 
from these principles are often called the Mind! in plate, or more properly the 
Reissner-Mindlin plate models. 

Problems with constraints need special treatment. The constrained variational 
form seeks a solution in a subset in which the constraint is identically 
satisfied. The discretized version of this form finds a solution in an 
approximation subset where the constraint is also identically satisfied. The 
construction of such a finite element basis is not always straightforward. A 
number of approximation procedures have been developed for plates and shells 
(often referred to as the discrete Kirchhoff method) and for incompressible 
problems. 

When inequality constraints such as contact and friction conditions are 
imposed on solid and structural problems, they are formulated by variational 
inequalities. The formulation is useful to demonstrate the mathematical 
complexity of this class of problems. Lions and Stampaccia (1967) present the 
derivation without directly referring to minimization theory. They also 
discuss the existence and uniqueness of a class of variational inequalities. A 
tutorial is prepared by Oden, Kikuchi (1977) and Kikuchi, Oden (1979), which 
includes a discussion on the finite element approximations of the inequality 
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constraints. For a general class of inelastic problems, the Hu-Washizu 
variational formulation can be cast into a system involving an inequality and 
equalities. This modification appeared in Nagtegaal and Nakazawa (1985). Again 
difficulties are encountered when an analyst attempts to construct the 
approximation subset directly for this class of problems. 

The mixed variational formulation and its finite element counterpart were 
first derived for problems with incompressibility constraints. Hermann (1965) 
introduces the Lagrange multiplier for the equality constraint and sets up 
finite element equations which include the displacement and discretized 
Lagrange multiplier as unknown variables. The unconstrained displacement 
subspace is used in this computation in conjunction with the subspace in which 
the multiplier lies. 

The perturbed form of the Lagrangian functional satisfies the constraint 
approximately. The Lagrange multiplier can be eliminated from this form to 
produce the penalty function form. Zienkiewicz (1974) provides documentation 
for the motivation and a preliminary discussion on the penalty form. The 
variational argument for the constant dilatation element in Nagteqaal, Parks 
and Rice (1974) also results in the same formulation for the problem of 
plasticity. As discussed in Oden, Kikuchi (1982) and Kikuchi, Song (1982), the 
same modification is possible for problems with inequality constraints. 

The elimination of the constraint in the set of admissible functions 
introduces a different problem related to the admissible space for the 
Lagrange multipliers, which looks equally if not more complicated. The 
difficulty is the condition for stability related to the space for the 
Lagrange multiplier. The continuous form satisfies the necessary condition for 
stability identically, but the finite element approximations do not satisfy 
the same condition and need special attention. This is discussed in detail in 
the next section. 

Identification of the Lagrange multiplier by the field variable leads to a 
simplified (i.e., modified) variational formulation originally derived by 
Rudiger (1960). Its numerical solution in the form of the hybrid displacement 
method is reported by Kikuchi and Ando (1972), which was the subject of 
controversy. A general procedure is discussed in Zienkiewicz (1977) and is 
applied to combine the solution of the boundary and finite element procedures 
[Zienkiewicz, Kelly and Bettess (1979)]. Zienkiewicz and Nakazawa (1984) 
demonstrate the utilization of such a principle to problems with an interface. 
Oden and Martins (1984) applied this procedure to generate a surface force 
term in variational formulation for friction problems. In recent papers, Manq, 
Gallagher (1977, 1978), Haugeneder, Mang (1979), and Mang, Hofstetter, 
Gallagher (1985) repeatedly discuss the virtues and vices of this simplified 
form. Some of these papers are philosophical rather than technical, and the 
statements made therein might be discounted. The conclusions drawn in these 
papers are negative, and do not provide a confident basis to fully exploit the 
simplified form for a wide range of applications. But this does not disqualify 
the concept. Xue and Atluri (1985) discuss the lack of stability of this class 
of variational formulations. However, the result does not fully explain the 
successes and failures reported in the previous papers. 
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Instead of directly referring to the variational functional, the finite 
element method can be derived from the method of weighted residuals [Finlayson 
(1972)]. Early literature on the method discusses the virtues and vices of 
these two formulations [Oden (1969A,B)]. In recent years textbooks have been 
published which develop the finite element concept almost exclusively from the 
Galerkin method concept [Fairweather (1978) and Fletcher (1984)]. The slight 
discrepancy between these two notions is sometimes useful (but not essential) 
when the method is applied to non-self-adjoint problems often encountered in 
fluid mechanics. Examples are use of the Petrov-Galerkin methods to generate 
the upwinding effect for convection dominated problems [Hughes and Brooks 
(1979) among many others]. 

As discussed in the standard mathematics textbooks [Lions, Magnes (1972) and 
Aubin (1972)3, these two different ways of telling the story (i.e., the 
variational functional method and the weighted residual method) are 
essentially slightly different manifestations of the theory of differential 
equations. The use of the variational method argument in which the 
differential equations are weighted and integrated before discretization is 
more useful for gaining insight into the mathematical structure of the problem 
and the finite element approximation process. In modern literature, the word 
" variational " is used almost universally, including the weak formulation for 
non-self adjoint problems. The word " Galerkin ", which still appears in the 
literature, often causes confusion in communicating the mathematical reality 
behind the numerical methods and is to be avoided. 

3. 1.2. 2 Element Technology 

Improvement of element performance by virtue of mixed and hybrid variational 
formulations is currently the subject of active research. The original 
motivation to use these approaches was to overcome the difficulty of 
constructing compatible displacement finite element approximations for fourth 
order problems such as for Kirchhoff plates and shells. The construction of 
(^-continuous elements is extremely complicated, and the resulting basis 
functions are prohibitively expensive to compute. Argyris, Fried and Sharpf 
(1968) describe one of the highest order polynomial interpolations defined on 
triangles and used in practical computations. An incomplete cubic triangle was 
developed and later documented in Zienkiewicz (1971) in which integration by 
the area coordinate system is introduced. This particular element satisfies 
the slope continuity only at the nodes, and the approximate displacement field 
is discontinuous over the element edge. To justify the use of a wide range of 
incompatible elements, the patch test was proposed by Irons (1975) and used 
extensively as a tool to validate new elements. Irons and Loikkanen (1983) 
document the recent development of this test. Taylor, Simo, Zienkiewicz and 
Chan (1986) attempt to develop the test for a wider range of element 
construction such as the mixed method. The relaxation of the slope continuity 
requirement by these formulations allows the use of simple C°-continuous 
elements. Various techniques such as the mixed and hybrid methods have been 
developed for this purpose. Application of mixed and hybrid approaches has 
been extended beyond shell structures to include solid continua. 
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The development of element technology is concentrated on a simple class. The 
four-node quadrilateral and eight-node hexahederal elements are subjects of 
recent research and development. Due to the complexity of higher order 
interpolation functions, few if any high performance elements have emerqed 
which use more than quadratic functions. Exceptions are work on nine-node 
Lagrangian shell elements which often exhibit superb numerical performance. 

An extreme simplification of the approximate solution method is to use a 
piecewise constant displacement field as was investigated by Kawai (1978, 1980 
and 1983). The method produces poor results for the displacement field, but 
the quality of the stress solution and the accuracy of the limit load 
predicted by this type of method is quite good. 

Modern element technology utilizes the flexibility of the mixed and hybrid 
variational formulations and does not directly refer to the displacement 
principle. This paradigm shift legalizes "tricks" referred to as the 
"variational crimes" extensively used in finite element computations to 
increase numerical performance. The differentiation between the mixed and 
hybrid methods has become a subtle issue. Quite often these concepts are used 
to derive the stiffness equations which contain entries at precisely the same 
locations as the conventional displacement method. This is discussed from an 
historical perspective by one of the founders of the methodologies [Pian (1978 
and 1983)]. A series of monographs compiled by finite element developers 
indeed focuses upon the development of mixed and hybrid methods including an 
elaborate version of the differentiation issue [Glowinski, Rodin, Zienkiewicz 
(1979) and Atluri, Gallagher, Zienkiewicz (1983)]. A survey has been compiled 
by Atiuri, Tong and Murakawa (1983). In addition, progress is reviewed in a 
conference proceeding compiled by Spilker and Reed (1985). The name 
"nonstandard" finite element methods was given to this class of elements by 
Brezzi (1979) with a detailed mathematical analysis for fourth order equations. 

For the purpose of clarity, the mixed and hybrid finite element methods are 
defined as: 

Mixed methods - Variational forms and finite element equations employ 

multiple unknown field variables which are independently 
interpolated. Some of these variables may be eliminated 
algebraically from the finite element equations. 

Hybrid methods - Constraints are introduced in the finite element forms by 
virtue of Lagrange multipliers defined over the 
inter-element boundaries. These additional terms are 
often eliminated from the algebraic equations resulting 
in a form similar to conventional displacement models. 

In recent publications [Pian and Chen (1982), Pian and 
Sumihara (1984)], the definition of hybrid methods has 
drastically changed. The surface integral terms are 
dropped from the formulation resulting in mixed stress 
finite elements with piecewise discontinuous 
interpolation for stress. Often the strain interpolation 
function is generated from the stress approximation under 
a linear isotropic elastic assumption. 
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The subtle issue related to this class of finite elements involves numerical 
stability, often referred to as the Babuska-Brezzi condition [Babuska (1973) 
and Brezzi (1974)]. This condition holds for the continuous problem if the 
solution exists. This fundamental question on mixed and hybrid finite elements 
was raised by Oden (1982) in the context of the penalty function formulation 
with selective-reduced integration for a linear, incompressible problem. The 
mathematical foundation of the stability issue comes from the existence 
theorem of variational forms often referred to as the Lax-Mi Igram-Babuska 
lemma which first appeared in Babuska and Aziz (1972). 

The stability argument, as applied to constrained problems, has been addressed 
in many papers. Babuska, Oden and Lee (1978) introduced the stability criteria 
for a mixed hybrid method for linear second-order equations. For 
incompressible problems, Fortin (1977) studied the mathematical structure of a 
simple mixed formulation, e.g., bilinear displacement and piecewise constant 
pressure. A rate of convergence is detected in which the Babuska-Brezzi 
condition plays a crucial role. The analysis has been extended to include 
several other finite element interpolation schemes in Fortin and Thomasset 
(1979). A lecture by Girault and Raviart (1979) covers most of the important 
abstract results related to the finite element analysis of incompressible 
problems. Bercovier (1978) looked at the perturbed form of the mixed finite 
element equations and studied the role of stability condition. Oden, Jacquette 
(1984), Oden (1982A), Kikuchi, Oden, Song (1984), and Oden, Kikuchi, and Song 
(1982) discussed the necessity of finite elements having to satisfy this 
condition. The authors state that their results discredit popular selective 
integration elements such as the four-node linear and the nine-node quadratic 
Lagrangian bases with selectively reduced integration. From numerical 
experiments, a certain class of stable elements recommended by mathematicians 
(e.g., quadratics with single point pressure sampling) was found to be so 
extremely inaccurate that these elements are not usable in any practical 
analysis. Carey and Krishnan (1982) investigated the stability of popular 
penalty finite elements and redefined the stability estimates. Experimental 
observations are included to support the mathematical argument. 

Zienkiewicz and Nakazawa (1982) demonstrate the possibility of relaxing this 
condition in practical computations. Similar argument can be found in Malkus 
and Olson (1982 and 1984). These papers strongly suggest the possibility of 
using popular selectively reduced integration schemes for penalty function 
finite elements. A proper postprocessing scheme would satisfy a posteriori the 
Babuska-Brezzi condition. Johnson and Pitkaranta (1982) look into the same 
problem and draw a positive conclusion. Kikuchi, Oden and Song (1984) propose 
a stable recovery scheme for unstable selectively reduced integration 
elements. Zienkiewicz, Taylor and Baynham (1983) introduce practical quadratic 
elements which satisfy the condition a priori . In the Stokes flow setting 
Bratianu, Atluri and Ying (1984) develop a hybrid finite element and discuss 
stability in the same context. This documentation insists strongly on the 
necessity of imposing the stability condition. 
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Le Tallac (1980) addresses this issue in the nonlinear elasticity setting, and 
proposes a patch of conventional elements which satisfies the stability 
condition. Glowinski and Le Tallac (1982) examined the numerical performance 
of popular selectively reduced integration elements. The incompressibility 
constraint in terms of deformation gradient for nonlinear kinematics was dealt 
with by the augmented Lagrangian form in these papers. 

The same stability issue plays a significant role in problems with inequality 
constraints. Without directly using the variational inequality, binary switch 
algorithms for contact and friction have been developed [Francavilla and 
Zienkiewicz (1975)] and implemented into commercially available codes. 
Variational inequality was introduced to the finite element community, and a 
systematic effort has been mounted to study the convergence characteristics of 
approximate solutions for constrained problems [Glowinski (1979)]. The mixed 
and penalty forms of variational inequalities provide the condition for the 
approximation of the surface forces. Kikuchi (1982A,B) discusses the stability 
of the penalized form of the variational inequality and generates the 
classical nodally lumped nonlinear spring for contact and friction situations. 
In addition, a postprocessing scheme for surface force recovery is developed 
therein. Oden (1983) develops an approximation scheme for this class of 
problems. Further development of the finite element approximation for 
penalized variational inequalities is reported in Campos, Oden and Kikuchi 
(1982) including the stability considerations and smoothing algorithms. 

A simple test for the stability of mixed elements similar to the classical 
patch test is required by engineers. There is an ongoing research effort, but 
a robust method has not yet been reported in the literature. 

Selective and reduced integration schemes have evolved from finite element 
computation of plasticity. The reduced integration for continuum elements was 
published by Naylor (1974). A heuristic argument for a constant dilatation 
element, which is a four-node quadrilateral with piecewise constant pressure 
(i.e., one point integration), is found in Nagtegaal, Parks and Rice (1974). 

The global mixed scheme, in which the additional mixed variables are 
interpolated continuously almost everywhere in the finite element region, has 
not been used extensively for obvious reasons. The original implementation of 
the mixed method with C° continuous approximations for Lagrange multipliers 
involves the factorization of large indefinite systems of equations [Hermann 
(1965 and 1967)]. The mixed approximation is also used extensively in fluid 
mechanics applications [Taylor, Hood (1973) and Kawahara (1978)]. As mentioned 
earlier, certain global mixed interpolation schemes can be used in conjunction 
with the equilibrium iteration algorithm. The structure of this method and a 
number of fully worked out examples are contained in the report [Todd, 
Cassenti, Nakazawa and Banerjee (1984)]. 

In modern literature, however, the use of piecewise discontinuous 
interpolation is more common than the use of globally continuous mixed 
interpolation. The advantage of the discontinuous mixed interpolation in 
conjunction with equation solution by a frontal solution algorithm is 
demonstrated in Nakazawa (1984). 
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An example of high performance continuum elements based on the hybrid stress 
argument has been presented by Pian and Sumihara (1984). Further development 
was reported in Pian and Tian (1985). A systematic approach to improve the 
accuracy of linear continuum elements has been presented by Belytschko and 
Bachrach (1985), with a discussion of the underlying variational structure 
discussed separately by Stolarski and Belytschko (1985). The limitation 
principle discussed in this paper is an extension of the concept introduced by 
Fraejis De Veubeke (1965). 

The utilization of the Hu-Washizu principle for a class of assumed strain 
finite elements for plates and shells is discussed by Simo and Hughes (1985). 
Their argument not only extends the possible mixed finite element forms for 
structures but also provides a rational basis for construction of stress 
recovery schemes. This may potentially improve the performance of modern 
elements used in structural analysis such as the one proposed by Hughes and 
Tezduyar (1981). 

It is interesting to note that the assumed strain elements have evolved from 
the research in selective integration techniques used for the Reissner-Mindlin 
plate model and the discrete Kirchhoff assumption for a class of isoparametric 
elements. The original papers were published as early as 1971 [Pawsey, Clough 
(1971) and Zienkiewicz, Taylor, Too (1971)]. A comprehensive survey has been 
compiled by Pugh, Hinton and Zienkiewicz (1978). 

The use of under integrated elements, either by selectively reduced or by 
uniformly reduced quadrature, has become a common exercise in large-scale 
applications. A major drawback of these elements is that the resultinq 
stiffness equations are often rank deficient. The number of zero energy modes 
is often larger than the number of physically meaningful rigid body modes. 
These additional zero energy modes are often identified as the hourglass modes 
for linear elements and the Esc her modes for certain quadratics [Hinton and 
Bicanic (1979)]. When selective integration is used for the solution of 
incompressible problems, similar symptoms appear in terms of the approximate 
pressure field which is often referred to as the checkerboard pattern. 

These numerical noises are mathematically identified as the extra entries in 
the kernel for the finite element stiffness equations. Flanagan and Belytschko 
(1981) identify the extra kernel vector for the four-noded linear element with 
one-point quadrature. This vector is then used to generate the stabilization 
matrix to suppress the hourglass mode. A simple scheme was developed to 
control the numerical noise for the selectively integrated shell elements 
[Belytschko, Tsay and Liu (1981)]. An extension of this procedure for the 
plate element with one-point quadrature is developed in Belytschko and Tsay 
(1983). 

The same numerical instability occurs in finite difference computations, and 
similar hourglass control schemes have also been devised for practical 
computations [Maenchen and Sack (1964)]. 
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It is insisted in Oden (1983) that the removal of the kernel a priori is not 
necessary since the spurious modes can be removed by a post-processing scheme. 
However, the factorization of a rank deficient system of equations is often 
impractical unless the stabilization is introduced. 

3.1.3 Solution Algorithms 

It has been the habit of finite element analysts to assemble the global 
stiffness equations and then to factorize them. The advantage of this direct 
solution strategy, as discussed in standard textbooks [Strang, Fix (1973) and 
Zienkiewicz (1977)]* is the ability to obtain nodal displacements for any 
structure unaffected by the conditioning of the matrix with a limited and 
predictable amount of arithmetic operations. An often useful feature of this 
approach is the fact that the factor of the stiffness equations is available 
for different loading conditions. 

For transient problems, implicit algorithms involving matrix solution have 
been the main vehicle for determining the nonlinear dynamic response of solids 
and structures. An overview is provided by Belytschko (1983), and stability 
behavior is discussed by Hughes (1983) in the same monograph. 

The solution strategy for the stiffness equations is to utilize the sparsity 
and the banded structure of the coefficient matrix. Little progress has been 
made in the past five years except for minor refinements and improvements of 
existing computational procedures such as band or skyline storage (profile) 
solvers and frontal solution methods. In particular, a complicated, yet 
efficient frontal process was obtained in combination with a skyline storage 
scheme by Thompson and Shimazaki (1980). The skyline storage strategy is 
documented by Taylor in Zienkiewicz (1977) and an unsymmetric version of the 
frontal solution by Hood (1976), both with we 11 -documented source programs. 
Also, in order to increase capacity, a partitioned storage scheme for the 
frontal solution has been proposed by Beer and Haas (1982). However, it is not 
yet clear if these types of schemes are appropriate for boosting performance 
in vector and/or parallel processor computers. 

Interest in the iterative solution of algebraic equations was revived to 
enhance the capacity as well as the speed of finite element solutions. The 
monograph compiled by Liu, Belytschko and Park (1984) represents the state of 
the art of the mid-1980s. Algorithms based on the preconditioned conjugate 
gradient method are most commonly used for systems of well-conditioned 
equations. The basic theory and its application to displacement and mixed 
finite element equations are fully documented by Axelsson (1976 and 1978). An 
effort to devise a simpler preconditioner, such as the element-by-element 
approximate factorization, has been reported by Hughes, Winget, Levit (1983), 
Muller, Hughes (1984), and Muller (1985). A similar technique which uses a 
simpler representation of the stiffness matrix was proposed by Rashid (1969), 
and limited success was reported for the solution of three-dimensional 
problems. A comprehensive survey of these classical works on iteration is 
provided by Elsawaf and Irons (1982). 
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A construction of iterative solution algorithms for symmetric, positive 
definite systems of equations is discussed in an early paper by Hestens and 
Stiefel (1952). The conjugate gradient and preconditioned conjugate gradient 
procedures are used extensively for solution of finite difference equations 
and optimization problems. 

For a class of linear problems with equality and inequality constraints, an 
iterative algorithm is proposed by Arrow, Hurwicz and Uzawa (1958) and 
referred to as Uzawa's algorithm. In this process, constraints are imposed by 
means of Lagrange multipliers. The algorithm has been extended independently 
by Hestens (1969) and Powell (1969). In comprehensive reports published by 
Rockafeller (1971 and 1974), the algorithm is called the augmented Lagrangian 
method. The basic characteristics of this scheme in a finite dimensional 
setting are discussed in detail by Hestens (1975). 

In the finite element literature, the same concept has been developed and 
utilized for the solution of constrained problems such as incompressible 
elasticity. Applications of initial strain method in an iterative fashion are 
reported by Argyris (1965) and Zienkiewicz, Valliappan (1969). The 
penalization process and its iterative improvement are discussed by Felippa 
(1978) with the imposition of nodal displacement constraints as a particular 
example. 

A comprehensive treatise on the augmented Lagrangian method in the context of 
finite element application has appeared [Fortin and Glowinski (1982)] in which 
recent developments are documented with a number of fully worked out examples. 

An example of the augmented solution with nonlinear equality constraints in 
nonlinear elasticity is demonstrated by Le Tallac (1980). In the field of 
optimization, the use of iterative solutions is a lot more common than direct 
factorization, and the experience gained there is transferable to finite 
element practice. An example of such an argument is found in Kikuchi (1979). 
The standard textbook [Luenberger (1984)] is a useful source for potentially 
efficient solution schemes. In addition, a family of algorithms based on 
dynamic and viscous relaxation has been investigated in the modern context 
[Felippa (1984) and Zienkiewicz, Loehner (1983)]. 

The method for tracing the equilibrium path of nonlinear deformation processes 
in an incremental iterative fashion was developed originally by Marcal and 
King (1967). An incremental process to deal with the problem of plasticity by 
generating the tangent stiffness over an infinitesimal load step was developed 
by Yamada, Yoshimura and Sakurai (1968); this process is not iterative and 
often leads to sizable accumulation of truncation error. 

The solution for incremental displacements is often constructed using 
Newton-Raphson iterations. This scheme requires the construction and 
factorization of the global Jacobian matrix for every displacement update. The 
modified Newton method does not need matrix assembly and solution but is a lot 
slower to converge. Methods based on quasi-Newton updates have started to 
appear in the finite element literature. Two important papers are by Matthies 
and Strang (1979) introducing the inverse BFGS update, and by Crisfield (1979) 
in which the secant Newton update is tested. 
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The use of the classical conjugate gradient algorithm for the solution of 
nonlinear solid mechanics problems was attempted with limited success by 
Biff le and Kreig (1981). 

Application of a quasi-Newton update for fluid mechanics computations has been 
discussed by Engelman, Strang and Bathe (1981) in the Navier-Stokes flow 
setting. The Broyden update is used for nonsymmetric systems of equations. For 
the solution of actual problems, however, the method is combined with the 
traditional Newton-Raphson scheme and the direct substitution algorithm. 

In these papers, it was not concluded which method is the best overall. As 
surveyed in Crisfield (1982), the success of schemes depends heavily on the 
characteristics of the specific problem being solved. 

In the context of mixed iterative solution procedures, the variable metric 
approach was found to be more powerful than the conventional Newton-Raphson 
update for the problem of plasticity [Nakazawa, Nagtegaal and Zienkiewicz 
(1985)]. In particular, the secant Newton implementation of the Davidon 
rank-one quasi-Newton method was found to be the most efficient. 

Automatic load adjustment schemes within the iterative displacement update 
have been developed to improve the performance of the nonlinear solution, in 
particular near bifurcation points in the equilibrium path. The algorithm 
originated by Riks (1979) was developed further by Crisfield (1980). A family 
of algorithms has also been developed to trace the complicated equilibrium 
path in an incremental iterative fashion [Bergan, Holand, Soreide (1979) and 
Bergan (1982)]. A survey is available [Crisfield (1982)] with a number of 
fully worked out algorithms and examples. An algorithm combining arc length 
and search procedures was developed by Crisfield (1983) with a limited number 
of examples. In addition, a series of numerical experiments has been carried 
out by Abdel Rhaman (1982) with detailed documentation of the algorithms. 

3.1.4 Postprocessing 

Numerical processes which operate on the nodal values of the finite element 
solution to produce information such as stresses and strains are regarded as 
postprocessing procedures. The definition of postprocessing and an attempt to 
unify the concepts involved were discussed for the first time in a systematic 
manner by Babuska and Miller (1984A,B,C). To obtain a smooth and potentially 
accurate stress field from the displacement finite element solution, a 
primitive form of postprocessing has been routinely performed. In fracture 
mechanics applications, the concept has also been used to evaluate the stress 
intensity factor and other quantities. 

As pointed out by Melosh (1963) in the early days of finite elements, the 
stress values obtained by direct substitution of nodal displacement values 
into the stress-displacement equation are much less accurate than the 
displacements. Indeed, as shown in Ciarlet and Raviart (1973), the maximum 
error in stress values in the displacement method is a full order greater than 
that of displacement in terms of the L 2 _ norm and the maximum norm. 
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Schemes to calculate accurate stress fields from the given displacement 
solution have been the subject of intensive research. Oden and Brauchili 
(1971) propose use of the conjugate approximation for stress which is used as 
the global filtering scheme for stress. As discussed in detail by Oden and 
Reddy (1976), the conjugate approximation could possibly improve the quality 
of the finite element stress field. Theoretically, the best approximation 
property is expected in the mean square sense. The scheme has not gained 
popularity due to the expense of generating another finite element solution 
for the conjugate basis function. An improved computational procedure was 
proposed by Oden and Reddy (1973). For problems with stress concentrations for 
example, the global Gram matrix used to generate the conjugate basis is 
partitioned and the computational effort is reduced. Using the least square 
form of the constraint equation, a family of filtering schemes can be 
constructed [Hinton (1968)]. A local-global smoothing scheme for derivatives 
is constructed by Hinton and Campbell (1974), and its performance for 
eight-node quadratics has been demonstrated in Hinton, Scott, Ricketts (1975) 
and Zienkiewicz, Hinton (1976). The textbook by Hinton and Owen (1977) 
documents the details of this type of numerical procedure. 

Motivated by the success of numerical filters developed for unstable pressure 
fields in the penalty finite element solution [Zienkiewicz, Nakazawa (1982) 
and Lee, Gresho, Sani (1979)], Owa (1982) investigated global stress smoothing 
schemes. The method, referred to as variational recovery , uses the weak 
variational form of the constraint equation and introduces the lumped mass 
approximation to generate the approximate Gram matrix. The scheme is 
economical since it does not involve a matrix inversion, and it is applicable 
to a wide range of problems without significant modification of the existing 
displacement finite element code [Nakazawa, Owa and Zienkiewicz (1984)]. The 
convergence of stress recovered by this variational process is optimal in the 
mean square sense as demonstrated by numerical experiments. 

An algorithm to generate nodal continuous stress fields has been developed and 
used in an iterative solution scheme for the mixed finite element equation by 
Cantin, Loubignac, Touzot (1978) and Loubignac, Cantin, Touzot (1977). The 
numerical performance of this scheme is not satisfactory. Cook (1982) improved 
the quality of the stress recovery algorithm by introducing the finite 
difference argument on patches of regular finite elements. 

Jeyachandrabose and Kinkhope (1984) propose a strain filtering scheme for the 
eight-node quadrilateral element. The implementation is based on a local 
filtering concept. The filtered strain is constructed explicitly in the matrix 
assembly process resulting in the so-called "assumed strain" element. In this 
type of algorithm, the improvement of strain convergence is only by a constant 
dictated by the limitation principle which states that the rate of convergence 
of mixed and hybrid methods in which the strain and stress are interpolated by 
element discontinuous functions is the same order as for the displacement 
method with the same displacement approximations. 
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3.1.5 Conclusions 

A substantial body of literature on the algorithmic and computational aspects 
of finite elements has appeared in the past several years. This survey covered 
major technical papers and reports published in the present decade and a few 
historically important contributions published earlier. The trend in research 
is to fine-tune the element technology and to streamline the solution 
procedure. The mixed and hybrid variational formulations play a central role 
in both of these aspects. These techniques are capable of handling 
large-scale, highly nonlinear problems in engineering because of recent 
developments in numerical methods and the significantly improved performance 
of modern hardware. 

Directions of future research and development, whose objective is to provide 
more accurate and efficient computational procedures, have become clearly 
visible from the literature surveyed herein. There have been, and will 
continue to be, two paths to take: 

1. For a given finite element mesh, extract more information more 
accurately, and 

2. For a given computer (constraint on the speed and capacity), 
accommodate larger meshes and crunch numbers in a feasible period of 
time. 

The spatial and temporal approximation strategy will allow realization of the 
first target. As discussed in this survey, the mixed and hybrid finite element 
concepts have been used to derive simple and highly accurate elements for two- 
and three-dimensional continua and three-dimensional shells. Further work to 
exploit the full potential of these "nonstandard" methodologies will improve 
the quality of the numerical solution. In particular, the utilization of 
continuous functions for the derivatives of field variables, i.e., 
displacement gradient, strain and stress, will potentially result in a robust 
numerical scheme. The stability and convergence rate (if the method at all 
converges) are yet to be determined. The mathematical structure of the general 
mixed finite element approximations should be an important subject of future 
research. In addition, the analysis of special cases such as singularities in 
strain and stress and material interfaces must be treated carefully in this 
type of methodology. Construction of approximation techniques must include 
these details. 

Reduction of core storage requirements and increases in speed partly rely on 
the element technology. The number of points where functions are evaluated is 
proportional to the number of arithmetic operations. Lower order integration 
and control of the kinematic modes associated with under- integration are 
crucial ingredients in the element technology. Further development and more 
numerical experimentation are required in this area. 


18 



A large portion of the computing time is consumed by the solution of linear 
algebraic equations in finite element analysis. An efficient solver is a key 
ingredient for overall efficiency. Iterative solvers should be investigated 
with particular applications to ill-conditioned problems such as those arisinq 
from the Timoshenko beam and Reissner-Mindlin shell theories. In nonlinear 
calculations, the Newton-Raphson process repeatedly assembles and factorizes 
the global stiffness equations which are the most expensive operations. Quasi- 
and secant-Newton processes will potentially decrease the cost of nonlinear 
solutions. The combination of an iterative linear equation solution together 
with nonlinear algorithms will be the subject of research and development in 
order to realize high performance software. The augmented Hu-Washizu 
formulation will provide a framework for the development of a family of 
preconditioned iteration algorithms. The key issue in this exercise is to 
construct a good preconditioner for the nonlinear iteration. In transient 
computations, the partitioned and operator splitting methods have similar 
characteristics, and investigations may result in an efficient finite element 
package. Comparison of high performance implicit schemes and explicit 
integrators should be made on a state-of-the-art coding basis. The recent 
investigations indicate that viscous and dynamic relaxation might result in an 
efficient solution package. This is also a subject of further evaluation. 

Postprocessing is a powerful method to increase the information extracted from 
a given finite element mesh. The augmented Hu-Washizu argument suggests that 
the result could be used iteratively to improve the overall quality of the 
finite element solution. The error estimator and indicator developed as a part 
of postprocessing are useful in avoiding unnecessary mesh-refinement, leading 
to an economical analysis strategy. The combination of an a posteriori error 
estimate and self-adaptive mesh refinement is useful for improvement of the 
numerical solution. The details of implementation will continue to be a 
subject of research and development. Comparisons of h- and p-version 
refinement are to be made in the nonlinear environment. The solution method 
for the enriched finite element model is not trivial. Here again, potential 
utilization of iterative solution algorithms should play a major role. The 
mixed finite element concept will provide here a wider range of possible 
tactics to refine the mesh than the straight displacement method. 

Modern computers exploit multiple processors in parallel or in vector form. 
Coding strategy to avoid excessive conditional jumping in loops will increase 
the throughput of these computers. In the equation solution phase, 
vectorizability and parallelizability are important points for evaluation. Not 
enough evidence has been accumulated to decide upon which solver is the most 
suitable for those machines. Currently, the machine dependency and code 
dependency of performance figures are unavoidable, and rationalization is 
required. 
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3.2 Global Solution Strategy and Element Technology 

3.2.1 Introduction 

Iterative methods are presented in this section for the solution of a class of 
indefinite systems of equations arising from the finite element discretization 
of mixed variational formulations. 

The main advantage of mixed methods is explicit integration of the constraint 
equations, which leads to stable and accurate numerical approximations in the 
finite element analysis. Disadvantages of the mixed approach are mainly due to 
the additional ( mixed ) variables resulting in a larger indefinite system of 
algebraic equations. Attempts have been made to solve the mixed finite element 
equations directly, which is, however, prohibitively expensive for practical 
computations. 

An iterative process will be constructed using the conventional displacement 
stiffness matrix (or its perturbed form with respect to the discretized 
constraint equations) as the preconditioner. A primitive version of this 
solution strategy was tested in Cantin, Loubignac and Touzot (1978) with 
application to problems of linear elasticity. The nodal continuous stress 
field was constructed based on the extrapolation process and was fed into the 
equilibrium iteration loop. The results were such that the accurate stress 
field was obtained, whereas the quality of the displacement solution 
deteriorated considerably. 

3.2.2 The Augmented Hu-Washizu Formulation 

Based upon the conventional virtual work statement in terms of displacement, a 
minimization problem of a functional is obtained: 


1 (U J 4/ u i,j °ijrs u r,s d ? 


which leads to a finite element equation: 

K u = F 


( 1 ) 


( 2 ) 


with K being the stiffness matrix. The total number of stress/strain sampling 
points in a given mesh is denoted by Ng. With application to inelastic 
analysis of continua, a vector of length (n^ + n) N g , where n is the 
number of dimensions, is required to store the state of stress and the 
deformation history. 
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The finite element method derived from the Hu-Washizu principle is augmented 
by the displacement method, equation (2), resulting in: 


K 0 8 

0 0 -C 


B 


T 





( 3 ) 


and this form is used to construct the iteration algorithms. 

The matrix B is the discrete gradient operator , the integration of which is 
carried out~in an element- by-element manner. The population is identical to 
the conventional strain-displacement matrix. The matrix C represents the 
strain (stress) projection operator which generates nodal strain (stress) 
values from the displacement gradient (or the stress sample D B u) calculated 
at certain sampling points. 

In the following algorithmic discussions, the possible interpolation functions 
for mixed variables are not limited to the element discontinuous polynomials. 
One of the advantages of the mixed method to be utilized in the present 
investigation is indeed utilization of continuous strain and stress 
approximations almost everywhere in the problem domain. 

3.2.3 Iterative Algorithms 

First, the implementation of Uzawa's algorithm is considered. Using the 
notation := to denote the assignment operator, the flowchart is: 

(1) Initialize the displacement and residual vector u: = 0 and R: = F 


(2) Displacement solution by factorizing the stiffness equations K u = R 


(3) Strain projection using the third equation in the Hu-Washizu e: = C~ T B T u 
finite element form (3) - - 


(4) Stress Recovery using the second equation in the Hu-Washizu s: = C” 1 D e 

finite element form (3) - - 

(5) Form the residual R: = F - B s 


(6) If convergent ( HR II less than the prescribed tolerance), then exit; else 
proceed to step (7) 
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(7) Displacement update by back substitution u: 


u + K" 1 R 


(8) Repeat from Step (3) 


In linear elastic calculations, the computational overhead to carry out this 
iteration is significant in terms of the storage requirements. The stress, 
strain and residual arrays have to be assigned explicitly, requiring nN^ + 

(n 2 + n)N s more words where Nh is the number of displacement nodes and 
N s is the number of stress nodes. In the above flowchart, it is assumed that 
the stiffness equations are factorized only once at Step (2) and that the 
factors reside in the memory until a new system of equations is assembled. The 
update in Step (7) only involves a back-substitution. The discrete gradient 
operator £ is formed as part of the residual recovery process, requiring 
additional arithmetic manipulations. 

The key ingredient for economy in the iterative algorithms is to diagonalize 
the discrete projection operator £. Then, a vector of length N s (instead of 
N § 2 ) needs to be added. The same array is repeatedly used 1/2 (n 2 + n) 
times, over the components of the strain vector. 

Provided that the particular mixed scheme gives a better solution than the 
preconditioner, the iterative strategy can be viewed as an inexpensive tool to 
improve the quality of the results obtained via the displacement finite 
element method. The iteration loop from Steps (3) to (8) is a postprocessing 
operation by virtue of the strain/stress projection. The operation can be 
truncated at any stage, and the residual calculated from the mixed stress 
approximation is then utilized as an error estimator of the finite element 
discretization. 

Line search , as incorporated in Uzawa's algorithm, seeks the new search 
direction by the orthogonality condition: 

a u R (u ♦ A a u) = 0 ( 4 ) 

where Au and R(u) are the iterative displacement update and the residual load 
vector associated with displacement ij respectively defined by: 

A u = K“* R ; R = F - B s (5) 

and displacement is updated by: 

u: = u + A a u (g) 

For linear problems, the parameter is directly calculated from a formula: 

a u T R (u) 

A (7) 

A U B A s 

where As is the iterative stress update associated with the displacement 
update Au. 


22 



The algorithm referred to 
schematically as: 

as the steepest descent method is described 


(1) 

Initial ize 

u: = 0 and R: = F 


(2) 

Displacement solution by factorizing the displacement stiffness 
equation 

Ku = R 

(3) 

Strain projection 

e: = C' T B u 


(4) 

Stress recovery 

s: = C‘ T D e 


(5) 

Form the residual 

R*: = F - B s 


(6) 

If convergent ( II R II 
proceed to step (7) 

less than prescribed tolerance), then exit; 

el se 

(7) 

Iterative displacement update by back substitution a u: a K"* R 


(8) 

Strain projection 

e: a C T B (u + a u) 


(9) 

Stress recovery 

s: a C- T D e 


(10) 

Form residual 

R: a F - B s 




a u T R* 


(11) 

Line search 

A - T “ 

a U (R - R*) 


(12) 

Displacement update 

u: a u + A a u 


(13) 

Repeat from step (2) 
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The coniuqate gradient is now introduced to search for the solution in a more 
efficient manner, in this algorithm, the directions are orthogonal with 
respect to the symmetric positive definite matrix appropriate for mixed finite 
elements. Note that the global mixed stiffness does not need to be formed 
explicitly. The residual vector calculated in the mixed manner is used to 
reflect the variational structure. The flowchart is: 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 

(7) 

( 8 ) 


Initialize 


u: = 0 s: = 0, s = 1 and R: = F 


Displacement solution 

Strain projection 

Stress recovery 

Form the residual 

If convergent ( II R || 
proceed to step (7) 

Preconditioning 
Displacement update 


K u = R ; a u*: = u ; = u T R 

e: = C‘ T B u 

s: = C' 1 0 e 
R*: = F - B s 

less than prescribed tolerance), then exit; else 

a u: = K _1 R* ; = R* T a u 

0 : = / « 0 ; 
au*:=au+0au*; 

*0 : ' *1 


(9) Line search to find the search distance s with respect to a u* 


(10) Displacement update u: = u + A a u* 

(11) Repeat from step (3) 
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It is possible to use modern numerical methods to update the preconditioner in 
an effort to improve the convergence of the iterative processes. Typically, 
the utilization of quasi-Newton algorithms is an attractive alternative to the 
class of schemes discussed above. 

The particular algorithm considered here is the inverse BFGS update based on 
the product formula: 


1 n 

^ * n a * -j Vj T i k - 1 n (i * vj w/> (s) 

j=n k=l 


in which the update vectors are calculated from the residual based on the 
mixed finite element formulation as follows: 


where 


= 


A “j 




L' 



* £ Zi 

a £ ?j-i 



(9) 


( 10 ) 


The scheme is implemented with and without line search. Note that the 
summation convention is not implied in the above formulae for the subscripts 
representing the iteration counter. 


In the computational procedure, the updated inverse replaces the inverse of 
the preconditioning matrix that appeared in Step (7) of the flowcharts for 
Uzawa's algorithm and the method of optimal descent. 

In the implementation, a fixed number of update vectors is stored, typically 
10, for both linear and nonlinear computations. 


A variable metric alternative to the quasi-Newton method is the secant Newton 
procedure . The Davidon rank-one quasi-Newton update in the secant form is 
utilized, that is: 



+ 



( 


a u 


n-1 



In 



In 



( 11 ) 
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with n being the iteration counter, resulting in an update formula: 


where 


1 in * A a is * 8 4 vi * c 4 ;s • 


•* :s - ?n 


( 12 ) 


(13) 


and 


C * 


(* v 1 * 1 S - 1 U S-1 ) T ?n 
( A Vl * 1 “S - ‘“J-l) 1 In 


(H) 


A = 1 - C ; B = -C . (15) 


This procedure does not require the storage of update vectors. The line search 
operation can also be incorporated into the numerical implementation. The 
interaction of this update process and the mixed iteration is limited only to 
the replacement of the displacement preconditioning by the above formula. 

It is interesting to note that, according to the author's knowledge, no 
attempt so far has been reported on utilizing the variable metric process for 
the iterative solution of linear finite element equations. 

3.2.4 Nonlinear Solution Processes 

Consider a quasi-static analysis of infinitesimal deformation processes. The 
finite element solution strategy is to subdivide the loading history into a 
number of successive load increments and repeat the Newton type iteration at 
each load level. 

Algorithms developed here are to embed the mixed finite element concept in the 
Newton iteration. The conventional strategy of implementing iterative solution 
algorithms is to code a nested inner loop inside of the existing nonlinear 
equation solution. This strategy may be economical if the core storage 
requirement is reduced by the decreased population of the preconditioner, and 
if convergence is achieved within a reasonable number of iterations. 
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In the incremental processes, the role of the constitutive equation is to 
calculate the state of stress for a given strain increment and a strain level 
accumulated at the beginning of the increment. The construction of a material 
tangent is a computational issue, which possibly improves the displacement 
preconditioner in the nonlinear iterations. 

The concept inherent in displacement methods which requires the constitutive 
equation to be identically satisfied everywhere in the domain is not 
necessarily compatible with the actual computational procedure when it is 
applied to nonlinear problems. In addition, this concept requires that the 
material tangent be calculated accurately everywhere in the finite element 
mesh. The mixed finite element concept based upon the Hu-Washizu principle 
describes the control structure of most finite element codes in a more 
straightforward and perhaps rational manner. 

The variational statement for the displacement method assumes the 
constitutive equation to be satisfied identically so that the strain 
and stress can be eliminated from the calculations. The material 
modulus needs to be given exactly at the state of stress and given 
deformation history to be able to make such elimination possible. 
However, inelastic material response is often given only as an 
implicit function of deformation history, state of stress and 
internal variables such as back stress. It is not feasible, even for 
a simple elastic-ideal ly-plastic material over a finite deformation 
increment, to derive an appropriate explicit stress expression in 
terms of strain and other internal variables. Thus, the constitutive 
relation which is to be used in the construction of the stiffness 
equations based on the displacement method concept can not really be 
exact. 

The computational procedure for the state-of-the-art plasticity model is to 
include a subprocess to recover the stress based on the radial return concept 
in the iteration loop. The subprocess involves two steps: 

1. Calculate the trial stress assuming the current strain increment is 
elastic. 

2. If the process is plastic, determine the radial return direction and 
intersection with the yield surface and evaluate the total stress 
accordingly and return; or else if the process is elastic, return. 
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The nonlinear solution process at a given load level F and increment AF is as 
follows: 


(1) Initialize arrays 


(2) Construct 


(3) Update au: = au + du 


(4) Strain recovery Ae: = C” 1 B T au 


(5) Stress recovery: Use the subprocess defined above and subtract the s 
associated with the end of the previous increment to generate an 
incremental stress array as associated with the current displacement 
increment 


(6) Residual calculation R: = F + aF - B (s + as) 


(7) If convergent ( II R II less than prescribed tolerance), then exit; else 
repeat from step (2) 


The iterative solution converges when the equilibrium condition is satisfied, 
regardless of the quality of the tangent array Kj. Step (2) contributes to 
the process only as displacement preconditioning in the general nonlinear 
environment, even when the displacement type approach is utilized. 

The preconditioned iteration algorithm developed and tested in the present 
status report is limited to linear and nonlinear material problems 
undergoing infinitesimal deformation. It is noted that for problems of 
finite deformation, the same iterative procedures would apply. Finite 
deformation kinematics related to geometry updating can be built into the 
stress, strain and residual recovery calculations. Such geometry updating 
effects will be incorporated into the MHOST code in conjunction with the 
large displacement option to be developed during the fourth year of 
contract work. In the framework of the displacement finite element method, 
application of a preconditioned iterative solution without explicit 
recalculation of the stiffness matrix has been demonstrated by Crisfield 
(1980, 1982). Examples discussed in these references include bifurcation 
due to geometrical nonlinearity among others. 


au, du, Ae, as and R 
Ky du * F - R 
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It is interesting to note that the preconditioned iteration algorithm is 
similar to the Newton-Raphson algorithm commonly used for standard 
nonlinear finite element computations except for the variational recovery 
of nodal strain. In the Newton-Raphson algorithm, as applied to finite 
deformation processes involving geometry updating, the approximate tangent 
stiffness matrix is generated with respect to the configuration qiven by 
the last displacement update. The displacement vector is updated using 
this tangent stiffness matrix. Then the strain and stress recovery 
calculation is performed in the usual manner in which finite deformation 
kinematics is taken into account. The iteration continues until the 
residual vector becomes almost equal to the external force vector in the 
deformed configuration. It can be strongly argued that the consistently 
updated tangent stiffness matrix improves the convergence characteristics 
of the iterative process. However, as the algebra indicates, the final 
solution obtained by Newton-Raphson- like methods is governed only by the 
equilibrium condition with respect to the approximate stress field, 
regardless of the stiffness matrix update in the nonlinear calculations. 

The advantage of the Newton-Raphson method in which the tangent stiffness 
matrix is assembled taking the geometry update into account is that 
sometimes it converges in fewer iterations than other Newton-Raphson- like 
methods because the quality of the approximated stiffness matrix is 
somewhat better than those updated by algebraic means in the BFGS and 
secant Newton methods. However, it is questionable if the additional 
expense of reformulating the stiffness and factorizing it is justifiable 
for large three-dimensional computations. 

In comparison with the Newton-Raphson method, as briefly discussed above, 
the preconditioned iteration algorithms described in this status report 
are generally more economical than the method involving repeated stiffness 
matrix assembly and factorization at each iteration. The change in the 
residual vector contains sufficient information to update algebraically 
the tangent array or its factors and such an approach is used in the 
quasi-Newton algorithms implemented in the MHOST code. 

In medium to large-scale finite element computations, the cost of stress 
recovery is often comparable in order of magnitude to the factorization of the 
tangent stiffness matrix. It is often possible to save on computing costs by 
choosing an appropriate basis for the stress interpolation so that N s is 
less than N g . 

As discussed in the previous section, the numerical performance of the 
nonlinear solution is improved by introducing several numerical ingredients. 
Indeed the methods of quasi- and secant-Newton were originally developed to 
solve nonlinear systems of equations. 

Major operations in the iterative process are the update procedure of K-j- 
used in Step (2) and the displacement update in Step (3). An alternative to 
the modified-, quasi- and secant-Newton processes, for nonlinear problems is 
the Newton-Raphson update, which involves assembling and factorizing the 
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tangent stiffness array at every iteration cycle. This method is not a 
particularly attractive option unless a identical to the mixed stiffness 
matrix could explicitly be calculated. The line search is also not an 
attractive option since multiple evaluations of the stress factor are required. 

Use of the arc-length method to control simultaneously the incremental load 
level and the iterative process is a useful option in conjunction with the 
modified-Newton, constant metric method. Exactly the same procedure as the 
displacement type is usable and is expressed schematically by the following 
flowchart: 

(1) Set the residual vector R: = 0, initialize the incremental displacement 
vector au: = 0 . If the first increment, initialize the total load 
factor X. 

(2) Project the displacement by 
du: = Kj 1 F 


and if it is the first cycle of iteration, calculate the arc length and 
initialize the incremental load factor a\ and dX. Then, calculate the 
correction vector. 

du: = K^ 1 (XF - R) 


(3) Update the displacement vector 


au: = au + du + dX du 


(4) Find the total load factor update 
X: s x + dX 

and the incremental load factor 


AX: = AX + dX 


based on the spherical path formulation. 

(5) Form the residual in the mixed manner and then check the convergence. If 
convergent, start the next increment, or else repeat from step (2). 
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3.2.5 Element Technology 

3.2.5. 1 Element Integration Procedures 

Linear and quadratic Lagrangian elements are included in the MHOST element 
library. A number of different integration schemes are implemented in the 
numerical solution procedures. The data base for the finite element model is 
constructed on a nodal basis with the exception of distributed loads which are 
stored in the element array. 

To achieve stability and coarse mesh accuracy in the mixed finite element 
solution, not all the terms in the equations are integrated at nodes. For the 
evaluation of the stiffness coefficients and the internal force vector, 
two-point Gaussian quadrature is used for linear elements and three-point 
quadrature for quadratics in each coordinate direction. A nodal quadrature 
(element discontinuous trapezoidal integration) is used for the strain 
recovery operations. All the stress components are evaluated at the global 
nodes. 

To extract near optimal numerical performance of the elements implemented in 
the MHOST code, the shear strain components are evaluated at the internal 
element centroid. 

All the integration procedures are internal operations designed to generate an 
accurate finite element solution and are invisible to the users as results are 
reported primarily at the nodal points. 

3. 2. 5. 2 Coordinate Transformations 

The coordinate transformation process used for the selective integration 
scheme has been modified to improve the performance of the bilinear Lagrangian 
elements used in MHOST. In particular, the polar decomposition for 
two-dimensional elements has been rewritten using the simple Cayley-Hamilton 
formula. 

For the two-dimensional Jacobian matrix J, the polar decomposition: 

J = (ROT) U (16) 


can be found from 


U 2 - I x ( U ) U + I 2 (U) I = 0 

where = U t U and I is the unit tensor. Noting that 

where U 2 = 1) and I is the unit tensor. Noting that 

(CH) = J T J = U 2 


(17) 


(18) 
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where = U" 1 " U 
and, hence. 


det (CH) * (det U 2 ) 


the following is obtained 


Thus, 


(CH) - I^U) U ♦ I 2 (U) I = 0 


U = (i^CH) ♦ 2 ^(CH))- 1/2 | CH ♦ VT 2 (CH) 1 ) 


( 19 ) 


( 20 ) 


( 21 ) 


This formula can be coded in an extremely simple manner. In addition, this 
approach is easily extendable to three-dimensional cases/ 

3.3 Local Analysis Procedures 

3.3.1 Introduction 

This section discusses numerical solution procedures specifically designed for 
structures with embedded discontinuities such as holes and cracks. Two methods 
have been proposed and analyzed in the framework of the mixed iterative finite 
element solution. Both methods utilize the enriched space of strain 
approximations in the discretized Hu-Washizu variational formulation. The 
first method refines the mesh locally within the conventional finite element 
context and extracts information by performing an independent solution. This 
method is referred to as the subelement scheme. The second method includes 
singular functions in the strain expansions around pre-assigned singular 
points such as a stationary crack tip. This method is called the special 
function approach. 

The subelement refinement strategy enhances the finite element approximations 
by subdividing the global elements used in the global analysis. The 
compatibility of the subelement displacement field with the global 
displacement field is maintained only at global nodal points in the 
local-global iteration. A series of numerical tests is conducted which 
includes a linear-stress patch test to validate the concept. Both linear and 
quadratic Lagrangian elements have been introduced for subelement 
computations. However, the validation effort is concentrated on linear 
subelement schemes including a number of problems involving large stress 
concentrations. A limited number of numerical experiments utilizing quadratic 
subelements has also been conducted. 
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In subelement regions, all spatially varying functions such as temperature and 
material properties are interpolated by the same shape functions as those used 
for the displacement field. When quadratics are used for subelement 
computations, the nodal values of the global element are interpolated 
quadratically. In the solution process here, the algorithm also generates 
quadratic fields internally for strain, stress and material moduli. 

The theoretical development of the special function approach is summarized in 
the next subsection. In the following, the form of the regular finite element 
shape functions is left unspecified, the derivation is general enough to allow 
incorporation of these results with all mixed element types implemented in the 
MHOST program. 

In the special function development, the singular functions are introduced 
only to represent strain and stress fields. For the region in which special 
functions are invoked, other variables such as coordinates, displacements, 
temperatures and material properties are interpolated by the global regular 
finite element shape functions. 

A discussion of the subelement method is provided in the previous Annual 
Report (NASA CR-175060). A description of the development of the special 
function finite element method follows. 

The dimension of the displacement approximation subspace remains unchanged, 
and the additional computational resources required for structural analysis 
with embedded discontinuities are kept insignificant in comparison with 
standard displacement procedures. Examples included in this report demonstrate 
the validity and numerical performance of the proposed procedures. 

3.3.2 Crack Tip Singularity and Special Function Expansion 

In the vicinity of a crack tip, the finite energy assumption leads to a 
hypothesis on the behavior of strain energy W, that is 


U = a., e.. — l -F( 0 )/r as r — —0 

I J I J 


( 22 ) 


under the infinitesimal deformation framework, with g and e denoting the 
stress and strain tensors respectively. In equation (22), function F(0) is a 
non-singular function depending on the behavior of the structure. 

For a linear elastic material, the linear relation of stress and strain 
implies that: 


ij 


Fi.(e) 


.-1/2 


(23A) 


and 


e 


ij 


V* 1 r 


- 1/2 


{ 23 B ) 
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The behavior of an elastic-perfectly plastic material can be written as: 


e ra — — He) r' 1 
re 


( 24 ) ! 


with the asymptotic behavior of individual stress components being potentially 
well above the yield level, whereas the asymptotic behavior of the Mises 
stress itself takes the yield value. 

When the strain-hardening of the elastic-plastic material is described by a 
power-law model, the stress and strain exhibit Hutchinson-Rice-Rosengren (HRR) 
singularities [Hutchinson (1968), Rice and Rosengren (1967)], i.e.. 


e 


ij 





and 


'ij 


G i j ( o) 


N 

TFT 


(25) I 


(26) ; 


where and Gjj are nonsingular tensor functions characterizing the 
effect of global deformation. The HRR form of the singularities contains the 
linear elastic and perfectly plastic material response as subclasses. 


Introducing singular behavior in the strain and stress fields, the following 
is obtained: 


and 


U 



* V 91 r ' p 


(26A) 



.. + F . .(e) r 


ij 


U 


(26B) 


where the superscripts R and S denote the regular and singular parts of the 
strain and stress respectively. Note that the condition: 

p + q ■ 1 (27) 

is necessary to satisfy the finite energy requirement. 
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Mow consider the displacement field. The fact that the derivatives of 
displacement contain a singularity does not imply that the function itself is 
singular. Therefore, the following local statement is valid: 


R + S 1 /. + .. 

e ij * Mj tj “ 7 “l ,j u j,i 


If the virtual stress is decomposed as: 


x , R + . S 

5a s So + 4a 


the global displacement gradient strain expression is constructed as: 

f50s • [(.•: - IT ( U • . + U. . ) J dx a 0, 

1 J 1 J ^ 1 » J J * ■ 

where 0 is the problem domain. 

Substitution of equations (26) and (29) into (30) results in: 

c / R R * * R S 4.5 R 4.0 S v j 

^ (5o ij e ij do ij e ij 6a ij e ij 6a ij e ij } 


/ 6ff ij (u i,J + U J.i } d - X 


1 C 

■y / 5a,*- (u. • + U. .•) dx 


1J 1,J J»1 


The finite element process has been derived for the weak form of the 
strain-displacement equation (30). 

The regular terms of the strain and stress tensors are expanded by the nodal 
shape functions M (same form as the displacement shape functions) as follows: 


R M R R M R 

a ij s \ a ijK * e ij = K E ijK ! 
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where the uppercase subscripts are used for the nodal indices. Introducing the 
nodeless variable set representing the intensity of the singularity, i.e., 
e s and£ S , the singular terms are approximated by: 


•ij ■ r ' q F u (9i * r 


q m k ( © ) 


i jK 


(33) 


’ij 


- r 


-P 


“ij* 81 * 


.-P 


M., ( 9 ) 


i jK 


Substituting equations (32) and (33) to (31), the following matrix equation is 
obtained: 




(34) 


where 


and 


C* R = / M T M dx 
^2 — — — 


C* = / (r' p ) M T M* dx 

~ RS n - - - 


c* a / (r" q ) M* T M dx 
« SR Q ~ ~ 


* , - 1 *T * 

C cc a / (r i ) M M dx 

.S5 n ~ 


B R = / M T ( V N ) dx 

n - - - 


B S = / (r' q ) M* T ( VN) dx 
a ~ 


(35) 

(36) 

(37) 

(38) 

(39) 

(40) 
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with N being the nodal interpolation function for the displacement. 


Except for the case of linear elasticity and small-scale yielding in the 
vicinity of the crack tip, p = q = 1/2 does not hold, and the resulting 
coefficient matrix (the singular strain projection operator) is nonsymmetric. 

To extract the information associated with the singular behavior of the 
structure, the lumping technique cannot be used, and the full projection 
operation is inverted at the element level. 

! The stress recovery process is now considered. Symbolically, the local 

I stress-strain law is written in an incremental form as follows: 

I 

a = 0 e ( 41 ) 

For a linear elastic material, this expression automatically yields the proper 
I singular stress form, equation (26), with: 

R n R • s „ -1/2 n r , Q v 

I a=u£ ,a=r D G(e) (42) 

I 


For the more general case of a nonlinear stress-strain relationship, equation 
(41) is used in an incremental manner. Although the order of singularity for 
stress and strain can differ, the incrementally linear approach will still 
produce the desired singularity effects. This is because the stiffness matrix 
D used to relate stress and strain increments is determined by previous 
deformation and stress histories and thus contains the proper effects to 
maintain the required stress and strain singularity orders. 

For a power law hardening material, equations (25) and (26), the general 
stress-strain relationship takes the form: 


a = , (43) 

Here, a strain singularity of order (1/N+l) in r automatically generates a 
singularity of order (N/N+l) in stress. It is also interesting to note that 
for a perfectly plastic material, N=0 and p=l, q=0; while for a linear elastic 
material, N=1 and p=q=l/2. 

The above observations indicate that the conventional incremental stress 
recovery routine can be used at every point in the finite element domain. The 
proper singular strain input produces the proper singular stress field without 
directly referring to the singular stress interpolation functions. 

The remaining task here is to actually construct the interpolation functions. 
The regular part of the finite element strain expansion is given by the linear 
Lagrangian basis, typically in two dimensions: 


H. . | (1 € > (1 + nj n) (44) 

I 
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where Sj and 17 ^ (i = 1 ,2,3,4) are the corner coordinates of a bi -square in 
the isoparametric plane. For the problems without singular points, the regular 
strain interpolations associated with the mixed finite element approach 
accurately approximate the field far better than the conventional displacement 
method. 

To introduce the singular terms into the finite element expansion of the 
strain field, the singular point is first preassigned at a corner of the 
quadrilateral element. The radial coordinate associated with a singular point 
assigned to the i th nodal point is represented in terms of the element 
isoparametric coordinates by: 

r 2 * (S - $.) 2 ♦ U - ni ) 2 , (45) 


while the tangential direction is given by: 

e * tan" * | ( n - n ^ ) / ( £ - ) j » ^ 

also in the isoparametric coordinate system. Then, the singular shape 
functions are defined as: 

M* = ( 8 A 2 ) (tt/ 2 - e) U/4 - e) 

= (4/v 2 ) e (ir/2 - e) (47) 

M-| = (8/it 2 J 9 (© - ir/4) 

for j=i+l, k=j+l, i=k+l and j=j-4 
if j>4, k*k-4 if k>4, 1=1-4 if !>4. 


3.3.3 Coding Strategy 

The singular strain function is introduced in the element strain calculation 
subprogram at the preassigned points as input data. No global smoothing will 
be carried out on the singular terms. The additional singular functions are 
generic and usable for both global and subelement solutions. 

The effect of the singular stress field can then be evaluated when the 
residual vector is recovered from the finite element stress solution. The 
stress recovery subprogram is executed for every element with preassigned 
singular points. 
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3.4 Computer Code Development and Validation 

3.4.1 Introduction 

In this section, the new options added to the MHOST finite element program 
during Task IVB are summarized. The number of algorithmic options has 
increased by including several alternative iterative solution methods and the 
optional frontal solution subsystem. The analysis options have been enhanced 
to include a composite laminate data reader and direct input capability for 
anisotropic elastic material properties. Internally, the control structure of 
the code, the common blocks and the finite element data base have been revised 
mainly to make the code more readable and debuggable. 

All the development work has been carried out on a PRIME 9955 at MARC Analysis 
Research Corporation. Procedures have been developed to create and test 
versions in the VAX/VMS, IBM/CMS and CRAY/COS environments. These procedures 
are executed when the installations at United Technologies Corporation/Pratt & 
Whitney and NASA-Lewis Research Center take place. The conversion takes 
typically one to two weeks including the testing. 

All the codes are compiled using FORTRAN-77 compilers. However, most of the 
source program is compatible with the previous FORTRAN standard. 

3.4.2 Solution Capabilities 

Version 3.2 of the MHOST program, which is the development version used in the 
Task IV computer program development effort, currently supports the following 
options and a limited number of linear and quadratic finite elements. All of 
the elements are usable in a fully mixed iterative manner, but the quadratic 
elements may only be employed in subelement regions. The conventional 
displacement solution algorithm is provided as a useful option with the 
default being the mixed iterative method. 

The programming effort is briefly summarized with regard to the control 
parameter keywords in alphabetical order: 

♦ANISOTROPY 

The anisotropic material property option is flagged by this parameter. The 
orientation vector for the material axis must be added in the model data 
section. The material properties along the material axis must be given by 
either the *DMATRIX option (for the continuum elements, type 3, 7, 10 and 11) 
or *LAMINATE option (for the shell element type 75). The anisotropic plastic 
response of a material is described by the user subroutine ANPLAS as 
documented in the MHOST Users' Manual. 

*BFGS 

The inverse BFGS update procedure is invoked by flagging this option parameter 
with the default iterative algorithm being the straightforward Newton-Raphson 
scheme. 
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♦BOUNDARY 


The nodal displacement constraints are imposed using a penalty approach. 
♦BUCKLE 

The initial stress matrix is formed, and an eigenvalue extraction is performed 
to obtain buckling modes. This option can be invoked at an arbitrary step of 
the incremental nonlinear solution process so as to detect the change in the 
buckling load due to inelastic response of the structure. 

♦COMPOSITE 

This parameter data line invokes the composite laminate analysis option for 
shell element type 75. This option resets the element parameters. The number 
of integration layers is reduced to one, and all the components of stress 
resultant are used as nodal variables. In the model data section, the user has 
to supply all the components of the nodal force-deformation matrix by invoking 
the laminate option. 

♦CONJUGATE-GRADIENT 

The conjugate gradient method is invoked by adding this parameter. The line 
search option is automatically turned on. Note that this option cannot be used 
with other iterative algorithms such as the BFGS and secant Newton-methods. 

♦CONSTITUTIVE 

Three different constitutive formulations are included in the code for 
describing material behavior. They involve secant elasticity (simplified 
plasticity) in which the material tangent is generated for use with 
Newton-Raphson type iterative algorithms, von Mises plasticity with the 
associated flow rule treated by using the radial return algorithm, and the 
nonlinear viscoplastic model developed by Walker, in which an initial stress 
iteration using the elastic stiffness is utilized. A linear elasticity option 
can be flagged for experimental purposes, and the default is the conventional 
von Mises plasticity model. Anisotropic plasticity is invoked by updating the 
user subroutine ANPLAS. The dummy subroutine supplied with the MHOST code is 
always required to produce correct results for isotropic cases. 

♦CREEP 

Creep effects are taken into account by integrating the time history in an 
explicit manner. An optional self-adaptive time step size control algorithm is 
also available. 
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*D X S PL ACEMENTMETHOD 

This option invokes the conventional displacement model in which the residual 
load is evaluated directly at the integration points. In the residual 
calculation procedure, the element strain is calculated directly from the 
nodal displacements via the element strain-displacement matrix. Then the 
element stress is obtained as a product of the element strain and the 
stress-strain matrix defined at nodes and interpolated to the integration 
points. This operation is consistent with the formulation of displacement 
stiffness equations and no residual force is generated in the equilibrium 
iteration. 

For the linear elastic stress analysis, no iteration would take place when 
this option is flagged. In inelastic analyses, the material tangent is 
interpolated and multiplied by the integration point strain which is directly 
sampled at the quadrature point. This option cannot be used for the advanced 
constitutive model since the correct material tangent is not generated. 

*DISTRIBUTELOAD 

Body force and surface traction loadings are referred to as distributed loads 
in the MHOST program. The body force option includes gravity acceleration 
definable in any direction and centrifugal loading with the centerline and 
angular velocity specified by the user. The detailed description is available 
in the MHOST User's Manual. 

♦DUPLICATENODE 

The continuity of stresses at nodal points can be broken by defining two nodal 
points at the same geometrical location and connecting them by this option 
which enforces compatibility of displacements only. This option is used to 
define the connections between generic modeling regions. 

♦DYNAMIC 

The generalized Newmark solution algorithm is entered by setting this flag. A 
simple adaptive time stepping algorithm is employed at the user's option. 

♦ELEMENTS 

The elements included in this version of the code are described in Table I. 
Core allocation is performed for the nodal and element-quantities on the basis 
of maximum storage space requirements among the types of elements specified in 
this option. All the element types must be specified here including those only 
appearing in the subelement regions. 

The element type, ITYPE, is chosen to avoid possible semantics conflicts with 
the MARC finite element code in which element type identification numbers 1 to 
98 have already been designated to various solid, structural and heat transfer 
elements. For instance, the four node shell element with 6 dof (degrees-of- 
freedom) per node based on the Reissner-Mindlin theory is referred to as 
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element type 75 in both MARC and MHOST. Also the two node Timoshenko beam 
element is referred to as type 98 in both codes. Nine-noded quadratic elements 
are rarely used in commercial codes, and thus these types of elements are 
designated as types 101, 102 and 103 in the MHOST code. Element type 
identification numbers over 100 are unlikely to be used in general purpose 
packages, and thus no conflicts are anticipated. 

The parameters defined in the table are used for core allocation purposes and 
do not necessarily correspond to physical reality. For instance seven words 
are reserved to define the nodal coordinate data of shell element type 75 by 
NELCRD. The first three words are used for the nodal point location in the 
global rectangular Cartesian coordinate system and the next three words are 
reserved for the components of the unit normal vector on the shell surface (at 
the node) generated internally. The seventh word is used to store the 
thickness of the shell at the node. 

The number of stress components stored for shell element type 75 specified by 
NELSTR consists of the generalized stress components preintegrated through the 
thickness including the in-plane twist component. 

The number of material property data entries NELCHR indicates how many values 
are needed in the *PR0PERTY model data block. Enough entries are allocated to 
store linear isotropic elastic constants, density and thermal expansion 
coefficient. To specify more complicated material response such as 
elastoplasticity and/or anisotropic elasticity, other parameter data lines 
must be added which invoke memory allocation for additional material data 
storage. 

The index for implementing the appropriate constitutive equation, JLAW, is 
used to tell the program exactly which stress and strain components need be 
considered as well as to select the form of the matrix which relates the 
stress and strain components. Other internal flags are used to specify the 
optional nonlinear constitutive models. 

Further discussion of these parameters, as well as the definition of other 
element characteristics, can be found in Section G of the MHOST User's Manual. 
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Table I 

MHOST Code, Version 3.2 Elements and Parameters 



Beam 

P. Stress 

P. Strain 

Axsym. 

Brick 

Shell 

I TYPE | 

98 

I 3 / 101 | 

11/ 102 | 

10/ 103 | 

7 

75 

NELCRD | 

3 

1 2 | 

2 1 

2 1 

3 

7 

NELNFR | 

3 

1 2 | 

2 1 

2 1 

3 

6 

NELNOD | 

2 

1 4 / 9 | 

4/9 | 

4/9 | 

8 

4 

NELSTR | 

1 

1 3 | 

4 1 

4 1 

6 

9 

! NELCHR | 

3 

1 5 | 

5 1 

5 1 

5 

5 

j NELINTI 

2 

1 4 / 9 | 

4/9 | 

4/9 | 

8 

4 

! NELLV | 

3 

1 3 | 

3 1 

3 1 

3 

4 

NELLAY | 

0 

1 1 1 

1 | 

1 | 

1 

5 

NDI | 

1 

1 2 | 

3 1 

3 1 

3 

2 

NSHEAR | 

0 

1 1 1 

1 | 

1 | 

3 

3 

JLAW | 

1 

1 2 | 

3 1 

4 1 

5 

6 


NELCRD Number of coordinate data per node. 

NELNFR Number of degrees-of -freedom per node. 

NELNOD Number of nodes per element. 

NELSTR Number of stress and strain components per node. 

NELCHR Number of material property data for the element. 

NELINT Number of 'full' integration points per element. 

NELLV Number of distributed load types per element. 

NELLAY Number of layers of integration through the thickness of 
the shell element. 

NDI Number of direct stress components. 

NSHEAR Number of shear stress components. 

JLAW Type of the constitutive equation. 


♦EMBED 

The subelement iteration technology is flagged by this option which signals 
the code to allocate the working storage for the subelement data in a 
hierarchical manner. The actual subelement mesh definition and the nodal- and 
element-data storage allocation take place when the individual subelements are 
defined. 

♦FORCES 

Concentrated nodal forces are defined and stored in a incremental manner. Core 
allocation takes place only when this option is invoked. 

♦FRONTALSOLUTION 

The frontal solution option for the quasi-static analysis is implemented in 
this version of the code. Out-of-core storage devices are utilized and, hence, 
capacity of the program is increased significantly. 
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*GMRS 


Generic modeling regions are defined as collections of elements that model 
geometrically parametrized parts of hot section components. Multiple generic 
modeling regions in a given mesh are connected using the duplicate node 
option. Different parameters are specified for each generic modeling region, 
and the input data can be prepared separately. Internally, the complex of the 
generic modeling regions is treated as a single mesh for the purpose of 
constructing and solving the finite element equations. A table is prepared to 
report results separately for each generic modeling region. 

♦LINE-SEARCH 

The line search option requires positive action by the user to turn it on. 

This algorithm may be used in conjunction with all the iteration algorithms 
available in the MHOST program. When the conjugate gradient iteration is used, 
the program automatically turns on the line search. 

♦LOUBIGNAC 

Parameters for the numerical quadrature used in the mixed iterative processes 
are defined in a very precise way. Full integration, selective integration, or 
selective integration with filtering can be chosen for construction of the 
stiffness matrix. For residual vector integration, full and reduced 
integration can be selected. The strain integration can be performed either by 
using uniformly reduced integration, trapezoidal integration with the reduced 
shear strain approximation or the previous quadrature with the filtering 
option. There are three parameters for this option to select integration 
schemes as summarized below: 


Parameter 1: 


Parameter 2: 


Parameter 3: 


1. Reduced integration for stress recovery 

2. Full integration for stress recovery 

3. (Default) Trapezoidal rule for stress recovery 

1. (Default) Full integration for residual force 
calculation 

2. Reduced integration for residual force calculation 

1. Full integration for stiffness matrix assembly 

2. Selective integration for stiffness matrix assembly 


3. (Default) Selective integration with element Cartesian 
coordinate transformation for stiffness matrix 
assembly. 

♦LOUBIGNAC 3 1 3 is the system default and recommended for most of the 
applications. 
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*MODAL 


The free vibration modes for linear elastic structures are extracted when this 
option is invoked. The subspace iteration technique is utilized, and a power 
shift option is included. 

*N0DES 

All the variables are defined and reported at nodal points. In the incremental 
processes, deformation and stress histories are integrated and stored only at 
the nodal points. Note that this architecture economizes storage substantially 
compared with fully integrated finite element displacement methods. 

*N0ECH0 

This option suppresses the echo print of the model data. 

♦OPTIMIZE 

A bandwidth optimization based on the Cuthill-McGee algorithm is applicable to 
in-core solution processes. No optimization is available for the frontal 
solver. 

♦PERIODICLOADING 

For transient calculations, nodal displacements and concentrated forces can be 
input as sinusoidal functions using this option. 

♦POST 

The postprocessing file which contains all the information supplied to and 
generated by the code are written to FORTRAN unit number 19. This file is 
formatted and can easily be manipulated by commercially available 
postprocessing packages with minor modifications. The header record of the 
file is designed to be compatible with the MARC post file which is processed 
by many finite element graphics packages. 

All the displacement and reaction force components written on the 
postprocessing file are in the global Cartesian coordinate system. When the 
user specifies nodal coordinate transformation, the displacement is printed in 
the user-specified coordinate system on the lineprinter output. Then the 
results are transformed internally to the global system before generating the 
postprocessing file. The postprocessing program package does not need to be 
told to recognize the user-specified coordinate transformation. 

♦PRINTSETS 

The report generation is carried out on a nodal point basis with element 
integration point options provided by interpolation using the shape function. 
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♦REPORT 


The frequency of the line printer report generation is now controlled by this 
option. The default is to print at every increment. 

♦RESTART 

This option tells MHOST to read the restart tape and set up the system to 
resume an analysis from the point where the tape is written. 

♦SCHEME 

Parameters that control the characteristics of the time integration operator 
are definable by the user. The default is the average acceleration algorithm 
commonly used in nonlinear dynamic finite element analysis. 

♦SECANT-NEWTON 

This parameter activates the secant-Newton implementation of the Davidon 
rank-one quasi-Newton algorithm. The core-storage requirement for this 
procedure is significantly less than that for the BFGS update, and execution 
is relatively fast. This option is recommended for inelastic analyses of 
sol id-continua. 

♦STRESS 

Boundary conditions for stress can be specified by the user as an option, 
although no mathematical justification is yet available for this type of 
constraint. Any stress component can be prescribed at any nodal point. Simple 
numerical tests have shown that inconsistent imposition of stress boundary 
conditions can lead to rapid divergence in the iterative process. 

♦TANGENT 

This option invokes the modif ied-Newton iteration procedure. The tangent 
matrix is updated until a user-specified iteration count greater than or equal 
to one is met. No updating occurs in subsequent iterations. The default 
iteration count is equal to one, and the procedure generated by this value is 
known as the Kyi method of modif ied-Newton iteration. This option has no 
effect when the BFGS process is employed. 

♦TEMPERATURE 

Nodal temperatures are read and used to generate thermal strains. These 
quantities are used also for the evaluation of creep strain and the 
integration of coupled creep-plasticity models such as Walker's model. 
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★THERMAL 


Temperature dependent material properties are evaluated when this option is 
invoked, and the appropriate user subroutine is provided to the system prior 
to execution. This operation is not necessary for the conventional creep model 
since temperature dependence can easily be incorporated into the user 
subroutine CRPLAW. 

★TRANSFORMATIONS 

Coordinate transformations at nodal points are specified by this option in 
which the angle of rotation is to be prepared by the user so that the code can 
generate the necessary transformation matrices. The postprocessing file does 
not support coordinate transformations at this time. 

★TYING 

Multiple degree-of-freedom constraint equations are specified by the user 
through this option. This option is more flexible than the duplicated node 
option since constraints may be applied to individual degrees-of-freedom. Note 
that the constraint equations are generated only for the displacement 
degrees-of-freedom. 

The following parameters are used to signal the MHOST proqran of the presence 
of specific user subroutines: 

★UBOUN *UHQ0K 

★UCOEF *UPRESS 

★UDERIV *UTEMP 

★UFORCE *UTHERM 

A brief description of each user subroutine is given below: 

UBOUN specifies the prescribed displacement as a function of time and/or 
increment. This subroutine is called once for every increment, and the loop 
over the total number of constrained degrees-of-freedom must be included. 

UHOQK specifies linear elastic material response varying with respect to node 
and temperature. This subroutine is called whenever the nodal stress is 
calculated. 

UCOEF specifies anisotropic thermal expansion coefficients which may be given 
as functions of temperature. This subroutine is called whenever the nodal 
initial strain is calculated. 

UPRESS is called at the beginning of every increment and specifies the 
distributed load varying with respect to tine and/or increment. 

UDERIV is an interface prepared for a user to add a special element in the 
form of the strain-displacement matrix. This subroutine is entered from the 
driver routine of the MHOST element library. 
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UTEMP specifies temperature dependent isotropic elasticity and thermal 
expansion. This subroutine is entered repeatedly at nodes from the linear 
elastic constitutive equation routine and the thermal strain routine. 

UFORCE is executed at the beginning of every increment and specifies the 
prescribed nodal force as a function of time and/or increment. The loop over 
the total number of points at which force vectors are given must be coded into 
this subroutine. 

UTHERM is called from the incremental analysis driver at the beginning of 
every increment to prescribe the nodal temperature field as a function of time 
and/or increment. The loop over the total number of nodes must be included. 

As summarized above, the analysis capabilities included in the MHOST code 
cover most of the needs for inelastic calculation of turbine engine hot 
section components. Presently, the code consists of around 40,000 lines of 
FORTRAN statements including extensive, self-explanatory comment lines in each 
subroutine. 

Table I summarizes the elements currently available in the MHOST code, Version 
3.2, and parameters used to internally define the element characteristics. 
Further details are available in the MHOST Users' Manual. 

3.4.3 Solution Algorithms 

The line search algorithm has been implemented into the MHOST code. This 
option is available with any iterative solution scheme implemented in the 
code. As shown in the following example, it is not clear if the combined 8FGS 
update - line search approach is the most efficient solution strategy in the 
present framework. 

As shown in Table II, the line search technique improved the convergence 
characteristics dramatically for the cantilever beam problem (Figure 1). 
Improvement of the solution by the line search was no longer obvious when the 
scheme was applied to a prototype problem for three-dimensional shell 
structures as shown in Table III and Figure 2. 

The inelastic response of the cantilever beam (Figure 1) is illustrated in 
Figure 3 and Table IV. The convergence behavior of the quasi-Newton, BFGS 
update was found to be superior to the Newton-Raphson type process, in which 
the tangent stiffness matrix was reformulated and factorized at every step of 
the iteration. It is important to note that the algorithm implemented here is 
far more efficient than Newton-Raphson type processes. 
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CM CO LO 


Table II 


Elastic Beam: Convergence with Respect to Iteration 


Relative Residual 


Constant Metric 


Variable Metric ( BFGS) 


Iteration 

Count 


0 

1 


Constant 

Update 

1.0 

0.97637 

0.50296 

0.25836 

0.13250 

0.06790 


Line 

Search 


1.0 

0.47754-6 

0.12064-6 

0.98808-8 

0.30286-8 

0.35292-9 


Constant 

Update 

1.0 

0.98028 

0.40637 

0.30279-2 

0.35640-4 

0.20975-7 


Line 

Search 


1.0 

0.11215 

0.22835-1 

0.13770-2 

0.27029-3 

0.16302-4 




r° 

1 L i i in iii 


1 L, 1 1 1 1111 


YOUNG'S MODULUS 1.0 + 7 
POISSON'S RATIO 0.3 


1000 

1000 


Figure 1 A Beam Problem - Linear Plastic Analysis 


49 



Table III 


A Three-Dimensional Shell Problem 
Convergence with Respect to Iteration 


Relative Residual 


Constant Metric Variable Metric 


Iteration 

Count 

Constant 

Update 

Line 

Search 

Conjugate 

Gradient 

BFGS 

Secant 

Newton 

0 

0.51222 

0.51222 

0.51222 

0.51222 

0.51222 


98.868 

98.868 

98.868 

98.868 

98.868 

1 

0.29682 

0.17956 

0.19521 

0.28914 

0.28322 


60.786 

38.035 

41.186 

60.811 

56.219 

2 

0.19051 

0.92278-1 

0.13825 

0.10905 

0.72094-1 


39.838 

19.521 

29.249 

23.421 

15.514 

3 

0.13140 

0.56205-1 

0.56362-1 

0.45218-1 

0.31623-1 


27.729 

11.990 

12.150 

9.7468 

6.8351 

4 

0.93775-1 

0.19583-1 

0.30214-1 

0.20974-1 

0.19135-1 


19.892 

4.1953 

6.4940 

4.5096 

4.1212 

5 

0.67875-1 

0.15912-1 

0.17466-1 

0.13668-1 

0.16991-1 


14.445 

3.4039 

3.7416 

2.9369 

3.6519 


A 



1 - TEST PROBLEM FOR ELEMENT TYPE 75 KINK ANGLES 

Figure 2 A Three-Dimensional Shell Problem 
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THE STRESS STRAIN LAW 
YOUNG'S MODULUS 1.0+6 
POISSON'S RATIO 0.3 
YIELD STRESS 7000. 

HARDENING PARAMETER 1.0 + 5 

THE LOAD DEFLECTION CURVE 

THE RESULT IS OBTAINED BY THE MODIFIED 
NEWTON - ARC LENGTH (SPHERICAL PATH) 
METHOD 



5.0 10.0 15.0 

TIP DEFLECTION 


Figure 3 A Beam Problem - Elastic Plastic Analysis 
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Table IV 


Solution of the Elastic Plastic Beam Problem 


Method 

Load 

Level 

Number of 
Iterations 

Normalized 
Tip Deflection 

Newton Raphson, Constant Update 

1.0 

5 

0.9917 

1.1 

2 

1.0997 


1.2 

10 

Fail to converge 

Newton Raphson, Line Search 

1.0 

1 

1.0000 

1.1 

1 

1.1000 


1.2 

3 

2.1937 


1.3 

10 

Fail to converge 

Modified Newton Line Search 

1.0 

1 

1.0000 


1.1 

1 

1.1000 


1.2 

9 

4.1340 


1.3 

3 

15.0842 


1.4 

5 

26.0958 


1.5 

5 

37.1663 

Conjugate Gradient 

1.0 

1 

1.0000 

1.1 

1 

1.1011 


1.2 

10 

Fail to converge 

Quasi-Newton/BFGS Update, 

1.0 

2 

1.0023 

Constant Update 

1.1 

1 

1.0985 

1.2 

7 

4.4550 


1.3 

5 

14.6930 


1.4 

9 

24.1375 


1.5 

5 

36.4625 

Quasi-Newton/BFGS Update, 

1.0 

1 

0.9997 

Line Search 

1.1 

1 

1.1016 


1.2 

3 

3.5301 


1.3 

10 

Fail to converge 

Secant Newton/Davidon Update, 

1.0 

2 

1.0002 

Constant Update 

1.1 

3 

1.1000 

1.2 

6 

4.4108 


1.3 

10 

18.5279 


1.4 

9 

24.7733 


1.5 

7 

34.7596 

Secant Newton/Davidon Update, 

1.0 

1 

1.0125 

Line Search 

1.1 

3 

1.1005 


1.2 

7 

5.8146 


1.3 

3 

16.7263 


1.4 

10 

Fail to converge 


Note: Parameters are: Maximum Number of Iterations 10, Relative Residual 

Error Tolerance 0.1, and Relative Displacement Error Tolerance 0.05. 
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The line search algorithm was then improved to reduce the number of residual 
load evaluations in the iterative process. In addition, the accuracy of the 
search distance parameter was improved significantly. The numerical results 
obtained by the new version of the code are compared with corresponding values 
from the previous scheme in Table V for the 20-element plane stress model of a 
cantilever beam. 


Table V 

Elastic Beam: Convergence of Modified Line Search, 
Conjugate Gradient and Secant Newton Procedures 


Iteration 

Line Search 

Line Search 

Conjugate 

Secant 

Count 

(Old) 

(New) 

Gradient 

Newton 

0 

1.0 

1.0 

1.0 

1.0 

1 

0.47754-6 

0.14709-7 

0.14709-7 

0.98772 

2 

0.12064-6 

0.51162-8 

0.51481-8 

0.16435-1 

3 

0.98808-8 

0.24806-8 

0.25119-8 

0.12846-3 

4 

0.30286-8 

0.12944-8 

0.13530-8 

0.38142-10 

5 

0.35292-9 

0.67459-9 

0.12350-9 

0.23028-10 


The preconditioned conjugate gradient algorithm has also been included in the 
MHOST code. The performance of this procedure was compared against the line 
search algorithm, as shown in the third column of Table V. Convergence was 
j extremely rapid with the residual monotonically decreasing. This procedure did 

i not, however, offer a particularly attractive alternative for solving 

inelastic problems. 

i The secant-Newton method of the form derived from the Davidon rank-one 

quasi-Newton update has also been implemented into the MHOST code. This option 
is particularly attractive when compared with the quasi-Newton method since 
storage of update vectors is not required. The performance of this approach is 
; included in Table I. The secant-Newton procedure converged faster than the 

i BFGS scheme for this application. 

I 

It should be noted that utilization of variable metric iterations in the 
linear equation solution for finite element analysis has not yet appeared in 
. the literature. 

| 3.4.4 Frontal Solution Subsystem 

Implementation of the frontal solution subsystem required a substantial coding 
effort that included a considerable amount of work associated with the 
reorganization of MHOST data structure. A core allocation routine was coded in 
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conjunction with a driver routine to control the quasi-static analysis by the 
frontal solution subsystem. The previously available band matrix solution has 
been left untouched and existing data decks can be used without any 
modification. 

The particular version of the frontal solver included in the MHOST code was 
designed for nonsymmetric indefinite systems of equations and utilizes 
diagonal pivoting within the front matrix. To avoid excessive read/write 
operations on a scratch file, a storage scheme was designed which utilizes the 
unused portion of work-space at the end of core-allocation as primary data 
storage. This space is referred to as the virtual disk file. Then the scratch 
file is assigned to secondary data storage, and read/write operations are 
carried out in a sequential manner. At the end of solution operations, the 
contents of the virtual disk are dumped to the scratch file, and this work 
space is released for other operations. This storage scheme is particularly 
efficient on machines without virtual memory but with high-speed secondary 
storage devices, such as the solid state device for CRAY computers. 

The frontal solution subsystem currently supports the following algorithmic 
options: quasi-static analysis, line search, conjugate gradient, tying, and 
duplicate nodes. 

A number of validation cases have been successfully executed on the PRIME 9950 
at MARC under the Primos 19.4 operating system, and the results are summarized 
in Table VI. It is noted that the code was compiled using the FORTRAN-77 
compiler prior to the testing, and it has also been compiled and tested with a 
FORTRAN-66 compiler. 

In Table VI, the memory requirement entries exclude the space used as virtual 
disk space. In the validation cases, the beam problem and the prototype vane 
example were suitable for band matrix solution since they consisted of 
regularly connected elements. Hence, the frontal solution shows a slight 
increase in computational effort, but the memory was less than the band matrix 
solution. The advantage of the frontal solution became significant when 
complicated meshes were processed, particularly for cases with tying 
constraints which increased the computational effort significantly for the 
band matrix solver. 

The program architecture related to the frontal solution subsystem is now 
summarized. The numerical treatment of duplicate nodes, tyings and 
transformations was significantly different from corresponding operations for 
these effects implemented in the band matrix subsystem. 
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Table VI 


Performance Tests for Frontal Solution Subsystem 




Memory 

Central Processing 

Unit Time 



Requirement 

(Seconds) 



Solver* 

(Words) 

Data Input 

Solution 

2-D Elements Model 

B 

10195 

0.952 

0.824 

For A Beam 

F 

9541 

0.942 

0.831 

Prototype Vane Model 

B 

46189 

1.100 

5.300 


F 

37559 

1.121 

5.670 

Plate With A Hole 

B 

12277 

1.000 

1.300 


F 

9617 

1.024 

1.103 

L-Shape Domain 

B 

27841 

1.645 

3.237 

(Continuous Stress) 

F 

18043 

1.712 

2.746 

L-Shape Domain 

B 

34459 (17404) 

1.782 

4.747 

(With Tying Option) 

F 

18695 ( 1626) 

1.839 

3.568 


I 

* B: Band Matrix 
i F: Frontal Solution 

The penalty function method was used to impose algebraic constraints of the 
duplicate node and tying types. Internally, an element consisting of two nodes 
was generated for each constraint data line in the connectivity data array. 

The stiffness coefficient was taken as 10^ times the value of the elastic 
modulus associated with the master node. Forward reduction was performed first 
for every regular element, and then these constraint elements were processed 
in a single stream of operations. 

The coordinate transformations were carried out on an element-by-element 
basis. When the element equations are assembled, a test is performed to 
determine if any one of the nodal points is subjected to the transformation 
and then the appropriate transformation operation is performed. 

For validation of these capabilities, a number of simple numerical tests were 
performed at MARC and the results were satisfactory, reproducing the same 
numbers as the band matrix option up to the fifth significant digit. 
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A MARC generated model of a three-dimensional burner blister specimen was also 
used to test the performance of the frontal solution subsystem. As shown in 
Figure 4, the model consists of 316 eight-noded hexahedral elements, 603 nodes 
and 123 nodal coordinate transformations. This example was also used to test 
the postprocessing file interface. The element renumbering to reduce the 
frontal matrix size was performed using an intuitive element renumbering 
scheme which involved geometrically sweeping the element identification from 
the center of the model to the outer circle. The statistics were: 

Band Solution Frontal Solution 


569,853 295,403 

157.8 201.5 


0.50719-3 

0.54622-3 


Core storage 
(32 bit-words) 

CPU time (seconds) 
per iteration 

Tip deflection 
(0 iteration) 

(4 iterations) 


The profile of the computed stress field has been successfully produced on the 
graphic postprocessing subsystem as shown in Figure 5. 



Figure 4 MARC Blister Specimen Model 
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BLISTER SPECIMEN 

2 28 64 97117 161 207 251 290 327 350 399 417 444 474 504 537 567 597 

NODES — • > « ■ ■ » » • » • » » • 



ARC LENGTH 


Figure 5 Mises Stress for Blister Specimen 
3.4.5 Anisotropic Material Property Reader 

At the request of the Contract Monitor at NASA-LeRC, a user-friendly format to 
deal with anisotropic material properties and composite laminate plates and 
shells has been added to the MHOST code. 

3.4.5. 1 Composite Laminate Option 

The composite laminate option has been completed, including a simple 
validation run. The material property data for the laminates are read directly 
into the elastic modulus array for shell elements in the format discussed in 
the paper Munich, Chamis (1975). For validation purposes, an example included 
in the above reference was solved. A cantilever beam made of a composite 
laminate is modeled using shell element type 75 in conjunction with *COMPOSITE 
and *LAMINATE options. 

The geometry and finite element model are illustrated in Figure 6 and the 
input data file is as follows: 

VALIDATION CASE FOR COMPOSITE LAMINATE OPTION 

C 

C EXAMPLE FOUND IN M.D.MINICH AND C.C. CHAMIS (1975) - 
C ANALYTICAL DISPLACEMENTS AND VIBRATIONS OF CANTILEVERED UNSYMMETRIC 
C FIBIRE COMPOSITE LAMINATES, NASA TECHNICAL MEMORANDUM NASA TMX- 

C 71699 

C 
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♦ELEMENT 8 
75 

♦COMPOSITE 
♦NODES 15 

♦BOUNDARY CONDITIONS 50 

♦FORCE 30 

♦END 

♦ITERATION 


5 

0.05 

0.00 

0.01 


♦COORDINATES 




1 

0.0 

-0.5 

0.0 

0.1 

2 

0.0 

0.0 

0.0 

0.1 

3 

0.0 

0.5 

0.0 

0.1 

4 

0.5 

-0.5 

0.0 

0.1 

5 

0.5 

0.0 

0.0 

0.1 

6 

0.5 

0.5 

0.0 

0.1 

7 

1.0 

-0.5 

0.0 

0.1 

8 

1.0 

0.0 

0.0 

0.1 

9 

1.0 

0.5 

0.0 

0.1 

10 

1.5 

-0.5 

0.0 

0.1 

11 

1.5 

0.0 

0.0 

0.1 

12 

1.5 

0.5 

0.0 

0.1 

13 

2.0 

-0.5 

0.0 

0.1 

14 

2.0 

0.0 

0.0 

0.1 

15 

2.0 

0.5 

0.0 

0.1 
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♦ELEMENTS 75 


1 

1 

4 

5 

2 

2 

2 

5 

6 

6 

3 

4 

7 

8 

5 

4 

5 

8 

9 

6 

5 

7 

10 

11 

8 

6 

8 

11 

12 

9 

7 

10 

13 

14 

11 

8 

11 

14 

15 

12 


♦PROPERTIES 75 

1 8 1.00000 0.00001 0.00001 0.00001 
♦LAMINATE 
C 

C CASE VI IN THE REFERENCE 
C 

1 15 


2315.E+3 

24.E+3 

O.E+3 

2315.E+3 

O.E+3 

82.E+3 

37.E+3 

37.E+3 

-55.E+3 

O.E+3 

O.E+3 

55.E+3 

O.E+3 

O.E+3 

O.E+3 

O.E+3 

1929. E+0 

20. E+0 

O.E+3 

1929. E+0 

O.E+O 

68. E+0 

31. E+0 

31. E+0 


♦BOUNDARYCONDITIONS 
1 1 

1 2 

1 3 

1 4 

1 5 
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1 6 

2 1 

2 2 

2 3 

2 4 

2 5 

2 6 

3 1 

3 2 

3 3 

3 4 

3 5 

3 6 

♦FORCE 

13 3 15.0 

♦PRINT 
TOTAL 
STRESS 
STRAIN 

♦END 

♦STOP 
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A 


Figure 6 Composite Beam Model Geometry and Finite Element Mesh 


The material is Thornel-75-S/epoxy with a fiber volume ratio of 0 . 6 . For the 
details of ply orientation, refer to the comment lines of the input data deck. 
Each nodal point has the same set of material properties as assigned by the 
*L AMI NATE data segment. After four iterations, a relative residual of 0.049835 
is reached and the tip displacements are: 



Oth Iteration 

4th Iteration 

Reference 

At Point A 

0.0888 

0.0929 

0.0887 

B 

0.0628 

0.0658 

0.0644 

C 

0.0421 

0.0444 

0.0434 


where Reference indicates the displacement results reported in Minich and 
Chamis (1975), who used the six-noded quadratic triangles in a 16-element 
45-node mesh. For comparison purposes, an uncoupled model was constructed, 

j- • i _ f _ 


resulting in deflection 

values of: 



2nd Iteration 

Reference 

At Point A 

0.0427 

0.0405 

B 

0.0211 

0.0212 

C 

0.0035 

0.0036 


61 



3. 4. 5. 2 Anisotropic Elastic Material Reader 

With the user subroutine UHOOK residing in the MHOST code, the user can 
specify any set of anisotropic elastic properties. An improved user interface 
has been developed which reads the entries of material modulus directly 
instead of invoking UHOOK. The additional model data parameter *DMATRIX 
invokes the subroutine DMATIN which interprets the numerical data discussed 
below. This option is available for continuum elements as well as for shell 
elements where the *C0MP0SITE and *LAMINATE options must also be invoked. 

The data format is: 

*DMATRIX This data line invokes the direct material modulus 

data reader at nodes. (NELSTR + 1) numeric data lines 
follow. 

The first data line includes two integers indicating the series of nodes at 
which the code will apply the material property data to follow: 

Integer 1: The first node in the series. 

Integer 2: The last node in the series. 

If the second integer is zero or left blank, the following material property 
data are applied only at the node specified by the first integer entry. 

Subsequent data lines consist of NELSTR real entries for each row of the D 
matrix at the nodes specified by the first data line. 

3.4.6 Coordinate Transformations 

The coordinate transformation process used for the selective integration 
scheme was modified to improve the performance of the linear Lagrangian 
elements used in MHOST. The polar decomposition formula for the 
two-dimensional Jacobian matrix was rewritten using the Cayley-Hamilton 
theorem. 

The numerical performance of the element implemented in the current version of 
the MHOST code was compared with that of the hybrid stress element proposed by 
Pian, Sumihara (1984), in which the element bending modes were dealt with in 
the isoparametric coordinate system. The example used here. Figure 7, was the 
pathological one discussed in this paper. Note that a=l means the angle of 
distortion is 45 degrees, the usual limit used in practical computations. For 
the range of distortion typically occurring in engineering computations, the 
preconditioning system did not perform as well as the rational hybrid stress 
element. However, the iterative algorithm generally outperformed the hybrid 
element as shown in Figure 7. 
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Figure 7 MHOST Selective Integration Versus Pian Rational Hybrid Stress 
3.4.7 Subelement Development 
Subelement Solution Procedure 


In the subelement region, the mixed finite element equations are constructed 
and solved iteratively. The displacement preconditioning is carried out in 
terms of the hierarchical nodal displacements in order to generate the 
incremental displacements in the subelement regions. Care is exercised to 
ensure that the effects of thermal and creep strains are correctly 
incorporated into the subelement load vector. 

Incremental strains are recovered at subelement nodes using a variational 
recovery process. Note that continuity of strain is no longer assumed over the 
global nodal points to which more than one subelement is connected. 

The constitutive equation is also integrated at subelement nodes using the 
subelement nodal strains. The contribution of thermal and creep strains is 
evaluated by using the subelement nodal temperatures interpolated from the 
global nodal temperatures specified by the user. The total and incremental 
stresses are stored at the subelement nodal points. 
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The subelement solution is coupled to the global finite element process 
through equilibrium iteration. The numerical integration procedure for the 
residual load vector plays an important role in controlling the overall 
behavior of the numerical solution. 

Two options are available to the user to control this coupling. The weak 
coupling scheme uses Gauss quadrature together with a simplified bilinear 
representation of the stress field. This scheme satisfies all the patch test 
requirements but its ability to extract local stress concentrations is not 
satisfactory. The strong coupling scheme, which does not satisfy the patch 
test for irregular meshes at this time, employs a special trapezoidal 
integration. The algorithm integrates the residual weighted by the area, which 
is conveniently available as the local nodal strain projection operator. 
Experience with the strong coupling scheme has led to modifications involving 
changes in weight factors. As the numerical example in the next section 
indicates, this modified strong coupling scheme detects stress concentrations 
effectively. However, the modified strong coupling scheme still slightly 
overestimates the residual when the contribution of the shear deformations of 
elements is significant. Further improvements are necessary and will be the 
subject of future investigations. 

The first validation problem consists of a corner singularity in a plane 
stress setting. Two global finite element meshes were generated. The coarse 
grid refers to the finite element mesh used in previous validation runs 
including 48 bilinear mixed Type 3 elements (Figure 8). The fine mesh was 
generated by subdividing all the elements in the coarse grid into 2 by 2 
bilinear elements. The subelement mesh refers to the coarse mesh with each of 
the three elements around the singular point represented by a 2 by 2 
subelement grid. 

At the singular point (re-entrant corner in the L-shape domain), the tying 
option is used to disconnect the strain and stress fields and, hence, the 
stress field can be approximated in an accurate fashion. 

The convergence characteristics are summarized in Table VII in terms of the 
relative residual. It is interesting to note that the BFGS update scheme 
diverges when applied to this problem. The solution at the singular point is 
summarized in Table VIII. A slight improvement of the nodal stress is observed 
in the subelement solution. 
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Figure 8 Problem Statement for Plane Stress Deformation of a Domain with 
a Singularity 
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Table VII 


Convergence of the Global Iteration for the 
Elastic Analysis of the L-Shape Domain 
(Prescribed Relative Residual Tolerance 0.01) 


Iteration 


Count 

Coarse Mesh 

Fine Mesh 

Subelement Mesh 

(I) Constant 

Metric Iteration 



0 

0.87591-1 

0.94754-1 

0.63856-1 

1 

0.48559-1 

0.48205-1 

0.41864-1 

2 

0.38770-1 

0.39204-1 

0.27499-1 

3 

0.33447-1 

0.34712-1 

0.18950-1 

4 

0.30079-1 

0.31803-1 

0.13658-1 

5 

0.27713-1 

0.29526-1 

0.10200-1 

6 

0.25930-1 

0.27699-1 

0.81305-2 

7 

0.24518-1 

0.26194-1 


8 

0.23360-1 

0.24518-1 


9 

0.22380-1 

0.23839-1 


10 

0.21531-1 

0.22893-1 


(II) Secant Newton (Davidon Rank One) 

Iteration 


0 

0.87591-1 


0.63856-1 

1 

0.48559-1 


0.41865-1 

2 

0.28219-1 


0.15284-1 

3 

0.26842-1 


0.11587-1 

4 

0.27726-1 


0.61043-2 

5 

0.24590-1 



6 

0.19135-1 



7 

0.18115-1 



8 

0.17086-1 



9 

0.14330-1 



10 

0.10645-1 
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Table VIII 


Stress and Displacement Solution: L-Shape Domain 


(I) Displacement at the Singular Point 


Coarse Mesh 
Fine Mesh 
Subelement Mesh 


X-Displacement 

-0.26425-3 

-0.25998-3 

-0.27157-3 


Y-Displacement 

-0.14922-3 

-0.13495-3 

-0.14972-3 


(II) Mises Equivalent Stress at the Singular Point. Note that duplicate nodes 
at the singular point result in different stress values at the same point. 


Global 

Global/Local 

Subelement 

(17484.) 

12539. 

11440. 

12589. „ 

B 

( ) refined mesh solution 


A 

C 


(11638.) 

(25634.) 

Global 

9296. 

19782. 

Global/local 

9316. 

19148. 

Subelement 

9621. 

20463. 


Table IX summarizes the results obtained by the modified integration scheme 
discussed above for both the global and subelement schemes. In the table, 
points A, B and C indicate the duplicate nodes at singular points in the 
bottom left, top left and bottom right elements, respectively. As indicated 
clearly, the stress values obtained by the subelement schemes are consistently 
higher than the finer mesh solution. This result qualitatively indicates that 
the modified integration scheme is now capable of detecting the stress field 
effectively and also properly transmits the subelement information to the 
global finite element solution. The quadratic subelement results show superior 
performance when compared to the linear subelement results. 

The above observations are preliminary and further quantitative numerical 
investigations need to be conducted. Possible code modifications to 
parameterize uniform subelement mesh generation will be considered to boost 
the performance of this procedure. 

The validation of the subelement code for three-dimensional analysis was 
carried out using a three-dimensional block under uniform mechanical and 
thermal loading as shown in Figure 9. The code produced the exact theoretical 
results for displacement, strain and stress. 
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Table IX 


Stresses at the Re-Entrant Corner 


Mises Equivalent Stress Value 



Point 

A 

Point 

B 

Point 

C 

Coarse Mesh 

9.296+3 

12.539+3 

19.782+3 

Fine Mesh 

11.638+3 

17.484+3 

25.634+3 

Coarse Mesh With Linear Subelement 
Global Node 
Subelement Node 

12.598+3 

14.085+3 

14.828+3 

18.468+3 

24.979+3 

30.901+3 

Coarse Mesh With Quadratic Subelement 
Global Node 
Subelement Node 

13.846+3 

14.949+3 

15.707+3 

22.391+3 

27.829+3 

35.750+3 




Figure 9 Three-Dimensional Subelement Validation 
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Quadratic Subelement Technology 

In this section, the validation of the quadratic subelements is discussed. To 
integrate the stiffness equations for quadratic Lagrangian elements, full 3 by 
3 Gaussian quadrature is implemented, whereas reduced 2 by 2 quadrature is 
used to recover the nodal strain. The validation cases reported here are 
limited to plane stress problems. In the range of numerical experiments 
discussed here, the above integration schemes are found to produce 
satisfactory numerical results. Further numerical study is necessary to 
construct an optimal integration procedure for quadratic elements. 

A series of two-dimensional numerical test cases was carried out for linear 
elastic and elastic-plastic problems. The results of a linear elastic beam 
analysis are summarized in Table X. The numerical model was the 20 plane 
stress element representation of the structure of which the two global 
elements on the fixed end were subdivided into both linear and quadratic 
subelements as shown in Figure 10. Note that exact values were reproduced at 
the nodes, although continuity of the displacement on the interelement 
boundary is not explicitly enforced except at the nodes. 

Table X 

Beam Problem: Subelement Validation 




Finite Element Solution 


Exact 


Linear 

Quadratic 

X-Coordinate 

Deflection 

Global 

Subelement 

Subelement 

0.0 

0.0 

0.0 

0.0 

0.0 

0.5 

0.0015 

(0.0060) 

(0.0030) 

0.0015 

1.0 

0.0060 

(0.0120) 

0.0060 

0.0060 

1.5 

0.0135 

(0.0180) 

(0.0150) 

0.0135 

2.0 

0.0240 

0.0240 

0.0240 

0.0240 


Note: Values indicated by ( ) are based on the nodal value calculated by the 
finite element interpolation functions. 



LINEAR SUBELEMENT 


QUADRATIC SUBELEMENT 


O GLOBAL NODAL POINT 
• SUBELEMENT NODAL POINT 

Figure 10 Subelement Refinement for the Beam Problem 
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The elastic deformation of a plate with a hole subjected to a uniform tensile 
load was studied using various finite element models. The basic mesh is shown 
in Figure 11 with the subelement division indicated by symbols. The stress 
distribution without any subelement refinement along the edge is shown in 
Figure 12 and the Mises equivalent strain contour plotted in Figure 13. 
Figures 14 and 15 represent the same information obtained from the uniformly 
linearly subdivided subelement mesh. The numerical values obtained from these 
analyses are summarized in Table XI. 


i 


HOLE IN PLATE UNDER IN PLANE LOAD - MHOST SETUP 
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12 3 4 5 6 

NODES □ □ O O Q O 

HOLE IN PLATE UNDER IN PLANE LOAD - MHOST SETUP 



O EQUIV. MISES TENSILE STRS 
13 1st COMP OF STRESS 
t 2nd COMP OF STRESS 

Figure 12 Stress Distribution Along the Edge A - B Global Mixed Solution 
with Coarse Mesh 



Figure 13 Mises Equivalent Stress Global Mixed Solution with Coarse Mesh 
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2 


3 


5 


6 


HOLE IN PLATE UNDER IN PLANE LOAD - MHOST SETUP 



O EQUIV MISES TENSILE STRS 
0 1st COMP OF STRESS 
^ 2nd COMP OF STRESS 


ARC LENGTH 


Stress Distribution Along the Edge A 
Linear-Quadratic Subelement Division 


- B: Coarse Mesh with 


1) 4.861+3 

2) 7.377 + 3 

3) 9.894 + 3 

4) 1.241 +4 

5) 1.492 + 4 

6) 1.744 + 4 

7) 1.996 + 4 

8) 2.247 + 4 

9) 2.499 + 4 

10) 2.750 + 4 

11) 3.002 + 4 



HOLE IN PLATE UNDER IN PLANE LOAD - HMOST SETUP 
EQUIV MISES TENSILE STRS 


Figure 15 


Mises Equivalent Stress: Coarse Mesh with Linear-Quadratic 
Subelement Division 



Table XI 


Elastic Stress Analysis of a Plate and a Hole: 
Stress and Strain at the Edge (Point A on Figure 10) 


Equivalent Mises Stress Horizontal Displacement 


Global Mixed Solution 


Coarse Mesh 

2.6026+4 

0.30068-3 

Fine Mesh 

2.9382+4 

0.34386-3 

Linear-Linear Subelement 

Global Node 

3.0054+4 

0.34912-3 

Subelement Node 

3.0059+4 

0.34912-3 

Quadratic-Linear Subelement 

Global Node 

3.0026+4 

0.35232-3 

Subelement Node 

3.0070+4 

0.35232-3 


Subelements for Embedded Holes 

A generalized version of the hole subelement has been fully implemented, and 
validation runs have been made. It is now possible to specify the number of 
subelement divisions in the radial direction (denoted by NUMRAD) and on the 
element edges (NUMEDG). Figure 16 shows the configuration of a subelement mesh 
for NUMRAD = 2 and NUMEDG = 2. The diameter of the hole and the subelenent 
type (either linear or quadratic) are the user-definable parameters for each 
hole subelement. 

The first validation case involved an analysis of stress concentration around 
a hole in a square plate. A uniform axial load was applied on the right hand 
edge, and the horizontal displacement was constrained on the left hand edqe. 
The vertical displacement was constrained at the left bottom corner. A uniform 
nine-element mesh was prepared as shown in Figure 17 in which the element at 
the center of the plate was subdivided into the subelement grid as shown in 
Figure 16. The result of this subelement solution was compared with the global 
finite element solution obtained from the reference mesh, as shown in Figure 
18. 
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Figure 18 Standard Finite Element Model with a Hole 


The initial solution of the model with an embedded hole was the uniform 
uniaxial state of stress with the displacement linearly increasing in the 
horizontal direction. In the residual force recovery phase of the mixed 
iterative operation, the subelements representing the hole produced 
out-of-balance forces at internal nodes. These forces represent the reduction 
in stiffness caused by the presence of the hole. The resultant forces were 
used in the equilibrium iteration to update the global nodal displacement. 

Stresses at points A and B of Figure 16 are presented in Table XII. The values 
are normalized by the uniform stress in a plate without a hole. As the 
tabulated values show, the subelement solution produced results consistent 
with corresponding results from the standard finite element solution. 
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Table XII 


Nodal Stress at the Edge of the Circular Hole 



Global Solution 

Subelement Solution 
Linear Quadratic 

Point A 


2.2615 

2.3198 

2.7087 

a y 

0.5916 

0.5977 

0.3929 

T xy 

0.0002 

0.0000 

0.0003 

Point B 

^x 

0.1752 

0.4339 

-0.2702 


-0.0813 

-0.2995 

-0.5958 

T xy 

0.0000 

0.0000 

0.0000 


Anisotropic Plasticity 

The generalization of the Mises yield criterion by Hill (1950) is implemented 
in the framework of the mixed, iterative solution approach. 

The anisotropic Mises equivalent stress with respect to the rectangular 
coordinates parallel to the material axes is defined by 

ai (°22 ' ° 33) 2 + a 2(°33 " a ll) 2 + a 3 ^l 1 * ° 22 ) 2 

+ 3a 4 o- 23 2 + 3a 5 <7 31 2 + 3a g <7.| 2 2 = 2a 2 


For plane stress problems, the formula is reduced to 


2 2 

a l°22 + a 2°’l 1 + a 3^ a l 1 


% 2 , 2 _ ,~2 
o"22' ^ 3a 6°12 ~ 2 ^ * 


where au (k = 1, 2 6) are defined as follows in terms of two vectors d 

and h which indicate the change in shape of the yield surface due to the 
anisotropic response of the material: 


a,= 


!_ . i i 

71*77 T? 
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do - 


1.1 i 

7 ? + 77 " ~1 

d 3 d i d 2 


and 


1.1 1 

71 + 7 ? ' 7 2 

d l d 2 d 3 


2 2 2 

a 4 = 7? a 5 = 71 a 6 = 71 
h 3 h 2 h l 


Note that isotropic Mises plasticity is recovered by setting each component of 
d = 1 and each component of h = 1 . 

The anisotropic formulation just described is suitable, for example, for 
modelling the behavior of sheet metal processed by rolling. This manufacturing 
process is likely to induce orthotropic anisotropy for the directions 
perpendicular to the roll. Further details on the physical assumptions 
involved in the generalized anisotropic plasticity model are presented in Hill 
(1950). Application of such a material model in finite element computations 
can be found in several commercially available general purpose finite element 
codes such as MARC. 


The cantilever beam problem, as defined in Figures 1 and 3 and loaded to 
undergo elastic-plastic deformation, is used as the validation case for this 
capability. An incremental load equal to 10% of the initial elastic load (i.e. 
the same incremental setting as used in the example summarized in Table IV) is 
applied for five consecutive increments and the secant-Newton/Davidon rank-one 
quasi-Newton update scheme is used in each increment. The material axes are 
assumed oriented parallel with the global coordinate directions. 

Under pure bending moment loading, the beam undergoes only elastic deformation 
up to load level 1.1. For the isotropic plasticity assumption, plastic 
behavior begins between load levels 1.1 and 1.2. Under these conditions, the 
x-direction stress component (oy) predominates at the locations on both the 
top and bottom surfaces of the beam, with the y-direction stress component 
(<7y) at these locations being smaller by one or two orders of magnitude. No 
shear stress occurs in this deformation process. To validate the 
implementation of the anisotropic plasticity option into the MHOST code, four 
sets of results have been generated. As a reference solution, the stress 
history for a standard isotropic plasticity model is calculated in a standard 
manner, the results of which are shown in Column 1 of Table XIII. In this 
table, the x-stress component (<r x ), as we ^ as the corresponding equivalent 
Mises stress (<f) values are given at load levels up through 1.5. The 
anisotropic plasticity option in the MHOST code is then exercised directly by 
making each component of both d and h take the value 1.0. As indicated by the 
results in Column 2 to Table XIII, isotropic plasticity is fully recovered 
using the anisotropic plasticity option with these special d and h vectors. 
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Material anisotropy is fully introduced in the remaining two sets of results. 
For Column 3 of Table XIII, the effective yield surface is extended by 50% in 
the global x-direction. This is accomplished by setting d-j = 1.5, so that 
the equivalent Mises stress values become much lower due to the reduced 
contribution of the predominant x-stress component ( ff x ). From Column 3 of 
Table XIII, it can be seen that the equivalent Mises stress value, even at the 
highest load level shown, i.e. 1.5, still has not reached the yield level. On 
the other hand, the results are far less sensitive to the effective yield 
surface being extended by 50% in the global y-direction. This is accomplished 
by setting d£ s 1.5 for the in-plane direction perpendicular to the beam. 

The extension by 50% in this direction alters the equivalent Mises stress 
values only very slightly as the load is increased. As can be seen from Column 
4 of Table XIII, the solution starts to deviate from the isotropic response at 
load level 1.4 but only ever so slightly in the fifth significant figure. 

Table XIII 

Isotropic and Anisotropic Plastic Response of a Cantilever Beam 

Anisotropic 


Load 

Level 


Isotropic 

Plasticity 

d = h = 1 

Plasticity 
d 1 = 1.5* 1 

d 2 * 1.5* 1 

1.0 

a x 

5998.8 

5998.8 

5998.8 

5998.8 


<7 

6002.9 

6002.9 

4001.9 

6002.9 

1.1 

°X 

6613.8 

6613.8 

6613.8 

6613.8 


a 

6568.2 

6568.2 

4379.3 

6568.2 

1.2 

°x 

7086.1 

7086.1 

7221.6 

7086.1 


a 

7001.4 

7001.4 

4768.4 

7001.4 

1.3 

a x 

7575.9 

7575.9 

7825.5 

7575.8 


a 

7556.9 

7556.9 

5162.8 

7556.9 

1.4 

°x 

7697.6 

7697.6 

8427.5 

7697.7 


a 

7601.2 

7601.2 

5559.8 

7601.1 

1.5 

°x 

7627.4 

7627.4 

9028.6 

7627.5 


0 

7673.1 

7673.1 

5958.2 

7673.1 


*1 Other entries of d and h are kept unity. 
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3.5 List of Symbols 


Alphabetical Symbol 
A, 8, C 
B 
C 

CH 

D 

D ijkl 

dx 

F 

F, 6 

h 

I 

I 

J 

K 

M 

M K 

N 

Ni 

n 

Ng 

Nd 

NS 


NOMENCLATURE 

Description 

Parameters in secant-Newton procedure 
Strain-displacement matrix 
Strain (or stress) projection operator 
Cayley-Hamilton matrix in polar decomposition 
Material modulus matrix 

Fourth order tensor component of material modulus 
Differential volume (under inteqral sign) 

Nodal force vector 
Function indicators 
Thickness for plates and shells 
Functional indicator 
Identity matrix (unit tensor) 

Jacobian matrix for isoparametric transformation 
Displacement stiffness matrix 
Element basis function vector 

Functional component of element basis function vector 

Plasticity power law hardening parameter 

Langrangian basis function associated with i-th node 
in an element 

Number of dimensions 

Total number of stress/strain samplinq points in a 
given mesh 

Number of displacement nodes 
Number of stress nodes 
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Alphabetical Symbol 
P 

P. q 
r 

R, R* 

ROT 

l 

s 

Jl 

Jl» Ji* 

Jij 

V 

£j 

u 

*j 

x 


Description 
Load at a given point 
Singularity powers 
Radial coordinate 
Residual force vector 
Rotation tensor 
Nodal stress vector 
Line search parameter 
Right stretch tensor 
Displacement vector 
Displacement vector component 
Space for admissible displacement variation 
Vector for BFGS update 
Strain energy density 
Vector for BFGS update 
Position vector 

Cartesian component of the position vector 


Greek Symbols 

& 

7 


5 0 » 5 1 
e 



Description 

Conjugate gradient parameter 

Coefficient parameter in plasticity power law form 

Vector differences in residual load vectors for BFGS 
approach ' 

Conjugate gradient parameters 

Strain tensor 

Strain tensor component 


80 


Greek Symbols 

6 

X 

v 

IT 

1 , 1 ? 

<7 

*10 

a 

da 

Sub- and Superscripts 
i » j » k, • • • 

1 y J) K 9 ••• 

MAX 

0 

S (Subscript) 

T (Subscript) 

R (Superscript) 

S (Superscript) 

T (Superscript) 

-1 (Superscript) 

Other 
( ), j 
( ’ ) 


Description 
Tangential coordinate 
Load factor in arc- length method 
Poisson's ratio 
Product indicator 
Isoparametric Coordinates 
Stress tensor 
Stress tensor component 
Problem domain 
Boundary of the domain 

Description 

Indices for vector and tensor component (when used 
as the subscripts); iteration and incrementation 
counter (when used as the superscripts) 

Nodal point counter 

Maximum number 

Initial quantities 

Shear and transverse shear component 
Tangent 

Regular part indication 
Singular part indication 
Indicates transpose of matrix 
Indicates inverse of matrix 


Derivative with respect to xj 
Derivative with respect to time 
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